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Abstract 


This paper investigates two mechanisms fundamental to sound 
generation in shocked flows: shock motion and shock deformation . 
Shock motion is modeled numerically by examining the interaction of a 
sound wave with a shock. The numerical approach is validated by 
comparison with results obtained by linear theory for a small- 
disturbance case. Analysis of the perturbation energy with Myers' 
energy corollary demonstrates that acoustic energy is generated by the 
interaction of acoustic disturbances with shocks . This analysis sug- 
gests that shock motion generates acoustic and entropy disturbance 
energy. Shock deformation is modeled numerically by examining the 
interaction of a vortex ring with a shock. These numerical simulations 
demonstrate the generation of both an acoustic wave and contact sur- 
faces. The acoustic wave spreads cylindrically. The sound intensity is 
highly directional and the sound pressure increases with increasing 
shock strength. The numerically determined relationship between the 
sound pressure and the Mach number is found to be consistent with 
experimental observations of shock noise. This consistency implies that 
a dominant physical process in the generation of shock noise is mod- 
eled in this study. 


Introduction 

This research was conducted to establish a better understanding of the physical nature of sound gen- 
eration in shocked flows. Increased understanding of the primary mechanisms in shock noise generation 
can be applied to a variety of aerodynamic problems. In helicopter and tilt-rotor aircraft, shock waves 
may form on blade surfaces and vortices from a preceding blade may interact with this shock, which 
results in impulsive noise. In supersonic engine inlets, combustion instabilities may result in shock 
oscillation that leads to what is commonly known as buzz (ref. 1). Subsonic transport aircraft often 
operate at conditions where the flow over the wing is transonic and unsteadiness in the flow over the 
wing can cause the shock to oscillate and generate noise. 

The primary focus during this research was sound generation in supersonic jets. When jet engines 
operate at supercritical nozzle pressure ratios, shock waves may form in the jet plume. Turbulence inter- 
acting with the shock waves generates high amplitude, broadband noise, typically known as broadband 
shock noise or shock-associated noise. The expected shock noise during the climb-to-cruise operating 
condition of a supersonic civil transport may ruin its chance of ever being put in production. Shock 
noise is an important design issue because of the effects on airport communities, aircraft interiors, and 
structural fatigue. 

Research toward the understanding of noise generation in jet flows has been ongoing for over 
40 years. Morley in an investigation of sound intensity in the far field of turbulent jets showed that the 
sound power is proportional to about the eighth power of the jet velocity (ref. 2). Lighthill in his pio- 
neering work on jet noise theory provided a theoretical basis for the eighth power law of Morley (ref. 3). 
Lighthill also provided a basis for an understanding of other jet noise phenomena, such as convective 
amplification (ref. 3). Lilley provides a review of classical jet noise theory and related experiments 
(ref. 4). There are two review papers specific to noise generation in supersonic jets that summarize sig- 
nificant contributions to the understanding of jet noise. The first, written from the perspective of an 
experimentalist, is the review paper by Seiner (ref. 5). The second, written from the perspective of a the- 
oretician, is the review paper by Tam (ref 6). The complicated nature of supersonic jet flow makes 
development of a comprehensive theory difficult and therefore the theories available for the prediction 
of shock noise in jets are largely empirical. 



The approach of research reported herein is to compute directly the sound generated by shock 
waves in supersonic jets. Advances in computer hardware that resulted in more computational speed, 
affordable memory, and parallel architectures combined with advances in software that resulted in more 
efficient, accurate algorithms, and better networking have made this approach feasible for relatively 
simple, two-dimensional problems. 

Direct computation of aerodynamically generated sound is useful because sound is inherently a 
component of a fluid flow field. The basic equations that govern sound are the same as those that govern 
fluid flow. For Newtonian fluid flows, these equations are the Navier-Stokes equations, or when viscous 
effects can be neglected, the Euler equations. As shown in figure 2, the acoustic portion of a fluid flow 
field consists of small perturbations on an inviscid, compressible flow. The generation of sound waves, 
however, often involves viscous nonlinear effects. Thus, an advantage of performing a direct simulation 
of the fully nonlinear equations is that both the sound generation and propagation are computed in the 
same analysis. In contrast, traditional acoustic methods require that the near field be found through a 
separate computation or experiment and the far field be computed with an acoustic analogy or Kirchhoff 
formulation. 
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Figure 2. Some equations of fluid mechanics. 
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Direct computation is currently not feasible for routine study of supersonic jet flow because the 
flow contains far too many scales to be resolved in a reasonable computation. Therefore, a combined 
modeling and direct computation approach is taken in this research. First, the complicated flow field of 
a supersonic jet is simplified into model problems, which are chosen to isolate mechanisms that are 
likely to generate sound. The modeling of supersonic jet flow is described in chapter 1. 

Chapter 1 begins with a description of features that are characteristic of supersonic jet noise and 
then proceeds to divide the complicated flow into several simpler problems that can readily be solved. 
Chapter 1 also presents an analytical model for shock noise generation. In this analysis, Lighthilfs 
equation is used to derive the monopole, dipole, and quadrupole source terms that are associated with a 
moving planar shock. The source terms are presented and then Green’s function is used to compute the 
far field sound that is associated with each of the source terms. All source terms were found to poten- 
tially contribute significantly to the far field sound. 

Chapter 2 addresses numerical methods and some challenges in the computation of flows with mov- 
ing shocks. Computation of sound generation in shocked flows is challenging because high-order accu- 
racy is required to determine the acoustic portion of the solution, while numerical dissipation is required 
at the shock to maintain stability. Most calculations presented in this report use a finite-volume imple- 
mentation of an essentially nonoscillatory (ENO) scheme. Because this algorithm is the basis for most 
calculations of this research, a brief description is included in this chapter. 

Once candidate sound generating mechanisms are identified and modeled and the numerical meth- 
ods are in place, the flow for each model problem can be computed and the acoustics extracted from the 
calculation. Chapter 3 describes the results obtained by an investigation of the first model problem: 
sound wave and shock wave interaction. Computations of the interaction of a sound wave with a shock 
wave in a quasi-one-dimensional nozzle show significant amplification of the sound wave pressure 
amplitude as a result of this interaction. The increase in amplitude of the acoustic pressure is shown to 
be a function of the Mach number upstream of the shock. Comparisons with linear theory are made for 
the small-disturbance calculations that validate the code. Results are also provided for the higher ampli- 
tude cases. In addition, an energy analysis is performed that shows that acoustic energy is generated 
during sound and shock interaction. 

Chapter 4 presents results that were obtained by the investigation of the interaction of a vortex ring 
with a shock wave. This sound generating mechanism is more complicated than the plane sound and 
shock interaction because the shock bends during the interaction, which generates alternating compres- 
sion and rarefaction regions along the acoustic wave front. Flow parameters downstream of the shock 
are observed and the conclusion is drawn that both acoustic waves and contact surfaces result from the 
interaction. Analysis of the results shows that the acoustic wave spreads cylindrically, that the sound 
intensity is highly directional, and that the sound pressure level increases significantly with increasing 
shock strength. The effect of shock strength on sound pressure level is consistent with experimental 
observations of shock noise, which indicates that the interaction of a ring vortex with a shock wave cor- 
rectly models some of the physics of shock noise generation. 

Figure Symbols 
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Chapter 1 

Modeling of Sound Generating Mechanisms in a Supersonic Jet 

The inherently complicated nature of jet flows makes both theoretical analysis and direct numerical 
simulation impractical for realistic flow Mach and Reynolds numbers. The structure of the flow field of 
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a supersonic jet consists of regions of laminar flow, turbulent flow, and transitional flow. Within the jet 
a myriad of structures of disparate scales exists, such as turbulent eddies and shock waves, that make 
analysis of the fluid dynamics practically impossible for general flows. Highly accurate numerical sim- 
ulation of these flows is impractical because of the disparity of the scales that must be resolved in the 
computation. 

To predict accurately the sound generated by complex jet flows, the essential elements of the jet 
fluid dynamics must be resolved in the computation. Experimental studies of noise generation in super- 
sonic jets have provided insight into what these essential elements may be. A result from an experimen- 
tal investigation of sound generation in a supersonic axisymmetric jet is included in figure 3. This 
photographic image was produced by a horizontal spark schlieren of an underexpanded supersonic jet 
and confirms the complex nature of the jet flow (ref. 5). As the flow exits the nozzle, it rapidly becomes 
turbulent. For the underexpanded case shown here, shock waves are also present in the flow and a Mach 
disk is evident approximately one jet diameter downstream of the nozzle. 



Figure 3. Schlieren photograph of underexpanded supersonic jet. 


For the acoustician, the wave field outside the jet plume is of particular interest. Several types of 
waves are present, which indicate that several sound generating mechanisms are responsible. Two types 
of wave fields are intense enough to be readily visualized on schlieren photographs. The first wave field 
appears to emanate from a region close to the nozzle exit near marker A on the figure and radiate at 
approximately ±30° from the jet axis. These waves propagate to the far field as acoustic waves. The 
waves are thought to be generated by supersonically convecting eddies that occur when the jet flow is 
heated and when the plume Mach number approaches two in unheated flows (ref. 5). The eddies create 
a form of sound known as eddy Mach wave radiation, which generally dominates the noise spectrum in 
directions of its dominant directivity. The primary directivity can be determined by computing the angle 
complementary to the Mach angle. The Mach angle is determined from the eddy convection velocity by 
}i = sin -1 ( MM) where M is the jet Mach number. 

The second type of wave outside the jet plume appears to emanate from the terminal locations of the 
shock waves in the mixing layer. One such wave field is shown emanating from a region in the proxim- 
ity of marker B in figure 3. This wave field appears to be more omnidirectional than the eddy Mach 
wave radiation because the wave fronts appear to spread spherically from the point of generation. This 
second type of acoustic wave is believed to be generated by the passage of turbulence through the shock 
waves and is referred to as broadband shock noise or shock noise. The broadband nature of this noise is 
illustrated in figure 4. 
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Figure 4. Spectrum of typical underexpanded supersonic jet. 


Figure 4 shows the spectrum of a typical underexpanded supersonic jet at 150° from the jet axis. 
The spectrum is characterized by three phenomena that are highlighted in the figure. The figure shows 
that these three phenomena tend to occur in distinct frequency bands. At low Strouhal numbers (a 
dimensionless frequency formed by the product of the frequency and the jet diameter divided by the 
fully expanded jet velocity), jet noise dominates. Jet noise is generated by turbulent mixing of the fluid 
within the jet. The peak in jet noise, at a Strouhal number of about 0. 1 , is believed to be generated by the 
large scale structures within the shear layer and is typically referred to as eddy Mach wave radiation. 
Both large and small scales within the jet generate mixing noise. The mixing of the fine scales produces 
the background level of jet noise. 

The peak of the spectrum is known as jet screech. Screech typically falls between the jet noise and 
the broadband shock noise in the spectrum. The accepted explanation for screech was originally pro- 
posed by Powell, who suggested that screech is generated when the turbulence interacting with shock 
waves in the flow develops a self-sustained aeroacoustic feedback loop (ref. 7). The resonant phenome- 
non, which generates a high amplitude sound at a particular frequency and its harmonics, is maintained 
by the shedding of a disturbance at the nozzle exit when sound passes the nozzle lip. Although Powell’s 
explanation of the screech phenomenon is useful in explaining general features of the phenomenon, 
there are features of jet screech that are still not understood (ref. 5). 

The high Strouhal number portion of the spectrum in figure 4 is called broadband shock noise. 
Broadband shock noise is present in all shocked jet flows and is due to the interaction of convecting dis- 
turbances with the shock waves in the jet plume. The broadband nature of this noise is due to the many 
scales of turbulent eddies in the flow. Experiments show that most of the broadband shock noise is 
directed slightly upstream (refs. 5 and 8). 

Although jet engines are designed to be shock free at design operating conditions, jet engines are 
often operated at off-design conditions. Thus, screech and broadband shock noise can contribute signif- 
icantly to the sound generated. 

Models 

Jet flows are too complicated to be practically analyzed or directly simulated. Therefore, the 
approach in this research is to identify elements essential to sound generation in supersonic jets and to 
analyze some models of these mechanisms for insight into the sound-generation processes. 
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Figure 5. Fundamental sound-generating mechanisms of supersonic jet flow. 


Experimental data indicate that noise is generated by different mechanisms within the jet. These 
sound-generating mechanisms associated with the complex flow of a supersonic jet are shown in 
figure 5. The first class of sound-generating mechanisms in the supersonic jet arises from the interaction 
of flow disturbances with the shock waves. The elements of the class that model the sound that is being 
generated by the interaction of flow disturbances with shock waves are illustrated in figure 5. These ele- 
ments are (1) the interaction of a vortex with a plane shock wave and (2) the interaction of a plane sound 
wave with a shock. These two elements will be discussed in detail in chapters 3 and 4. 

Additional models of sound generating mechanisms that are not directly related to the presence of 
shock waves are also illustrated in figure 5. The model presented for a wavy-plume structure is flow 
past a wavy wall. This models the component of sound generated by large scale structures within the jet 
plume. Clearly, Mach waves are expected to be generated by flow past a wavy wall, and in the model, 
these Mach waves are analogous to the eddy Mach wave radiation from supersonic jets. The vortex and 
vortex interaction model presented in figure 5 illustrates the interaction of turbulent structures within 
the jet, which are believed to be responsible for the background noise (jet noise) of supersonic jet flows. 


Evaluation of Sound Generation by Shock Motion With Lighthilfs Equation 

This section presents an analysis of sound generation by shock motion in the context of Lighthilfs 
acoustic analogy (ref. 3). 


Background 

Lighthill combined the equations that govern the conservation of mass and momentum into an 
inhomogeneous wave equation. A brief derivation of Lighthilfs equation follows. 
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The physical law governing the conservation of mass written in indicial form is 


dp , _ 0 

dt dxj ~ * 


j = 1,2,3 


( 1 ) 


where p is the fluid density, uj is the y'th component of velocity, t is time, Xj is the coordinate in the yth 
direction, and Q is the mass per unit time per unit volume injected into the fluid. 

The conservation of momentum equations are written as 


duj dP u 

p 37 + p “'5^ = -a^ + F ' iJ=1A3 


( 2 ) 


where F,- is an externally applied force, the stress tensor Py = p8y - Gy, where p is pressure, 6,y is the 
Kroneker delta, and Gy is the viscous stress tensor, Gy-2\ley — Xe kk 8y, where ey =1/2 (duj/dxj + 
dujfdxj), and |Li and X are the dynamic viscosity and second coefficient of viscosity, respectively. 

Taking the inner product of u i with equation (1) and adding it to equation (2) produces: 


3p«, d(pUjUj) dPjj 
dt dxj + dxj 


F i + U i<2 


(3) 


Differentiating equation (1) with time, subtracting the divergence of equation (3), and adding and 
subtracting c^V 2 p, produces Lighthill’s equation: 


^ p _ c 2 v 2 d = ag , tin 

dt 1 “ dt dx t dx i dx j 


(4) 


where is the free stream speed of sound and the stress tensor is defined as 


T ij= P U i U j - °ij + P 5 ij - C ~P S y ( 5 ) 

The terms on the right side of equation (4) are the acoustic source terms. The terms dQ/dt represent 
sound generated by unsteady mass addition; 3(F ( - + Q)/dx i represents unsteady forces; p u-uj represents 
nonlinear viscous effects, turbulence, and nonlinear propagation; and ( p - pc^)5.^. represents nonisen- 
tropic effects such as shock waves and heat addition. Effects of viscosity on the sound generation are 
represented by Gy. 


Source Terms 

Here, LighthilTs equation will be used to analyze the sound generated by shock oscillation. To sim- 
plify the analysis, dissipation effects are neglected and only the velocity component normal to the shock 
u is considered. Since there are no applied forces or unsteady mass addition, the source term simplifies 
to 


2 2 

T ii = p u + p-pc^ (6) 

The source term of Lighthill’s equation requires that the second partial derivative of T u be taken 
with respect to x. It is beneficial to use generalized functions because ordinary derivatives do not exist 
across the shock. The definition of a generalized derivative in one dimension is (ref. 9) 


d /(*,0 _ d /(*,0 

dx 3jc 


+ Afb(x-x s ) 


( 7 ) 
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where d/dx represents the generalized derivative operator, d/dx is the ordinary derivative, and A/ rep- 
resents the jump in the function /at the discontinuity located at x = x s . 

2 2 

Applying equation (7) twice to obtain d T u /dx , 


dT n 

_ ^ u 

dx 

dx 

?T n 

_ d 2 Tj 

dx 2 

dx 2 


+ at u &(x-x s ) 
3 


r + ^ X [AT i i 8(x - + A ("97) 5(jf _ x s) 


( 8 ) 


Thus, the source term in Lighthill’s equation for sound generated by a moving shock has compo- 
nents that resemble shearing stresses (quadrupole), unsteady forces (dipole), and mass addition (mono- 
pole), which are represented by first, second, and third terms on the right side of equation (8), 
respectively. 

In classical acoustics, where sound generation is considered in an ambient medium, the sound pres- 
sure associated with monopoles, dipoles, and quadrupoles is proportional to the second, third, and 
fourth power of the source Mach number, respectively. Therefore, the monopole source dominates at 
low Mach numbers, and the quadrupole source dominates at high Mach numbers. However, for sound 
generation at transonic speeds, each component may contribute significantly to the sound generation. It 
is a purpose of this section to analytically determine which term dominates the shock noise generation. 

To determine the behavior of the monopole, dipole, and quadrupole terms as functions of Mach 
number, consider a planar shock that is set in motion by an upstream disturbance in Mach number 
(fig. 6). The Mach number upstream of the shock is 

Afj = M q - M s cos((Ot - kx) (9) 

where Af 0 is the undisturbed upstream Mach number, M s is the maximum amplitude of the disturbance 
Mach number, co is the frequency of the disturbance, and k is the wave number. Disturbances in pres- 
sure, density, and sound speed upstream of the shock are neglected because these quantities are propor- 
tional to the square of the upstream disturbance Mach number. 


Monopole source term . The monopole term on the right hand side of equation (8) is 



( 10 ) 


y 



Figure 6. Shock oscillation used in determination of source term. 
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Considering for now only A(37"| j/Sjc), 



” 3jc( P 2 C 2 M 2 + ^2“P2 C ooj“0 r ^Pl <: l^l + P\ ~P\ C 

d(c 2 M 2 ) dp 2 


, 2 2 2 Y*P2 

~ | c 2 M 2 - c 


dM. 


,_ + 2 p 2C2M2 ___ + _^2p 1 c 1 M 11 - 

because pj , /?j, and C] are, to first order, constant. 


( 11 ) 


To evaluate equation (11) in terms of the upstream density, sound speed, and Mach number, the 
well-known shock jump relations are used. (See ref. 10.) These relations are presented for complete- 
ness. The second equality refers to an ideal gas with the ratio of specific heats, y = 1.4. 


P2 

2 Y m5-(y-1) 

7M^ - 1 

(12) 

Pi 

y + 1 

6 

P2 

(y + 1 )Mj 

6M^ 

(13) 

Pi 

(Y - 1 )Mj + 2 

M| + 5 

«2 

_ (y - 1 )M j + 2 

M( + 5 

(14) 

u l 

(Y+1)M^ 

6Mj 


To compute the monopole source term of Lighthilfs equation, the derivatives of the downstream 
pressure, density, and velocity are required. Taking the first derivative with respect to x , for the case of 
y = 1.4, results in 


dp 2 

dx 


7 dM j 5 

= 3 P ' M '~dT = 3 PlC|Ml 


3M, 

dx 


d(c 2 M 2 ) Cj k{M\ -5)3Af 1 
= ~6 


dx 

dp2 

dx 


M 


= 60p,- 


M, 


^5 + MjJ 


dx 


dM , 

2lT 


where = -M s & sin(co/ - kx). 

Substituting equations (12), (13), and (14) into equation (11) results in 



= -60 p j c 


2 

(«? * 5 ) 2 


sin(co/ - kx) 


(15) 

(16) 


(17) 


(18) 


11 



( 19 ) 


Substituting - M 0 - M s cos(o>/ - kx) and keeping only terms of order M s , results in 
Af-^ — J = bOpjC^ -sin(cof-fc;t) 


(^o + 5 ) 


Equation (19), normalized by 60p sin(co t-kx), is plotted as a function of M 0 (for constant 
wave number) in figure 6. 



Figure 7. Normalized monopole component of source term sound generated by sinusoidally oscillating shock. 


Two cases are shown in this figure. One curve represents the source term when the ratio of perturba- 
tion Mach number to upstream Mach number M/A/q is held fixed. The other curve represents the source 
term for a fixed perturbation Mach number M s . Clearly, the monopole source term of Lighthill’s 
equation reaches a maximum value. To determine the exact location at which it maximizes, the first 
derivative of equation (18) with respect to Mq is determined and set equal to zero. The value of M 0 for 
which this derivative is zero, when M is fixed, is found to be Af 0 = 75/3 = 1.29. The value of Afj for 
which (0/3MQ)A(9r u /9jc) = 0, when M/A/ 0 is fixed, is found to be M 0 = 75 = 2.24. 

Thus, the monopole source term reaches a maximum when, for constant shock velocity amplitude 
and constant disturbance wave number, A/ () = 75 / 3 . When the ratio of the shock velocity amplitude to 
the presoak Mach number is held constant, the source term maximizes at M 0 = J5 for a fixed wave 
number. 

For an acoustic disturbance, the wave number varies with upstream Mach number, 


k = 


(0 

cja+Afj) 

co cos (kx - co f) 


c,(l+M 0 ) ( 1 +M 0 ) 


M+o(MV 


2 a 


( 20 ) 
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The wave number may be approximated by k = co/cjO + Af 0 ) and substituted into equation (19) 
while maintaining a truncation error of the order of (M s ) 2 . Furthermore, if the sound speed approaching 
infinity is taken to be the stagnation sound speed, the relationship between c j and c ^ is 



Thus, the monopole source term for k = is 



-60p j c^coA/qA/^. sin(cor - /cjc) 



Equation (22) normalized by -60p\c oo <&M s sin(cot - kx) is plotted in figure 8. 


( 22 ) 
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Figure 8. Normalized strength of monopole component of source term for sinusoidally oscillating shock. 

Dipole source term . Consider now the dipole term of equation (8). This term can be written 

l[AT n b(x-x x )] = AT u S’(x-x s ) (23) 

since the jump across the shock is not a function of x. Thus, the jump in T\ j is of interest. 

AT n = p 2 M2 + P 2 -p2 c ~-(pl M l + / , I-Pl c i) ( 24 ) 

But, from the unsteady jump relations (see appendix A for derivation) 

p 2 (u 2 -u s ) = Pj(«j -u s ) 

P 2 U 2 ( U 2 - U S ) + P 2 = p 1 “ 1 ( M 1 + P\ ( 25 ) 

The change in pressure across the shock can be found with 

P 2 ~P\ = Pi(«i -« s )(«i -« 2 ) ( 26 ) 


13 



Substituting equation (26) into equation (24) and neglecting higher order terms in u s , results in 

~ -(P 2 - Pj) c oo (27) 

Therefore, upon substitution of equation (13) into equation (27) and simplification, the following 
expression for the dipole source term is obtained: 


AT,, = -5p jC^ 


M] + 5 


(28) 


Substituting Mj = Mq — M s cos(octf - kx ) into equation (28) and neglecting terms of order M and 
higher, results in 


ATn = 


- 5 Pl c o 


Mq + 5 L 


("J- 


1 )- 


12M s Mq cos(o)/ - kx) 


M 2 0 + 5 


(29) 


Unlike the monopole term, the dipole term contains both zeroth and first order terms in M s . Also, 
the wave number affects both the amplitude and the phase of the monopole term, while it affects only 
the phase of the dipole term. However, the monopole term and the first order component of the dipole 
term are identical except for the phase and the multiplicative factor of the wave number. 

In figure 9, there are no local extrema for the dipole source term when the flow is supersonic. This 
lack of extrema can be verified analytically by taking the derivative of equation (28) with respect to M 0 
as for the monopole case. 



Figure 9. Normalized strength of dipole component of source term for sinusoidally oscillating shock. 


Quadrupole source term. Finally, consider the quadrupole component of equation (8). The quadru- 
ple term requires that the second ordinary derivative of the stress tensor component T x l with respect to 
x be evaluated. Thus, 


?T U 

dx 2 


a, 

dJ 


^ 3w 

2p "S + “ 


2 3p dp 
3* + dx 



) 


_ (Tu 3 m dp fdu 2 2 d^p d z p 


(30) 
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The evaluation of equation (30) upstream of the shock is simple, since the derivatives in p, p, and c 
are high order and hence can be neglected. The evaluation of equation (30) upstream is 


,, d M \ 

~ 2p,C {l dx ) +W > dx 2 


Substituting = M 0 - M s cos(co/ - kx) and neglecting higher order terms in M s , this expression 
reduces to 

2 

dT,, -y 7 Af n Af cos(tDt -kx) 

j- = \0 Pi clk 2 -^ (32) 

dx i 5 + Af 2 

For an acoustic upstream disturbance for which the wave number is a function of Mach number, this 
becomes 


F u 2 MqM s cos(ou - kx) 

— = 2p,to 

c l f5 + M^(l+M 0 ) 2 


The evaluation of the quadrupole term downstream of the shock is more complicated since the sec- 
ond derivatives with respect to x of p 2 , w 2 , and p 2 are required. The second derivatives of the variables 
in equations ( 1 2), ( 1 3) , and ( 1 4) are 


5 2 p 2 60pj T 2 /dM i \ 2 / 3\0 2 A/ 1 

17 = r^ 5 - 3 "'(-ar) + ( 5M ' + "'br 

ax 5 + M, L 


P 2 5 2 ^ M \ 


d 2 “2 _ C 1 f ln r 3M lf (* M w 3 A 


where dM^/dx = -kM s sin((o/ - kx), and d 2 M 1 /dx 2 = k 2 M s cos(o )t - kx). Substituting equations (35), (34), 
and (36) into equation (30) provides 


2 f/3M n 2 d Af, 


= 2piCi — + M , — 

Hl 1 l dx 1 3 2 


12pic 


2 m\ + 20m\-25 /dAf,>2 


5 + Af 


3 \\ dx 


+ 1 2 p | c 


2 -5 Mj a Afj 

o° 9 9 

(5 + Af?l dx 


Substituting Af , = A/ 0 - M s cos(co/ - kx) and neglecting higher order terms results in 


f 1 1 2 2 (2 

■y 1 = 10p,cV 2_1_ m 2 -1 Jco8(fi)f-fcc) 
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Substituting in equation (20) for k = k(M 0 ), the quadrupole term becomes 

2 M 0 M s 

= lOpjto — (A^-Ocostcor-ih:) (39) 

2 (5 + Mq J 1 + Mq 

Note that the quadrupole source term contains only first order terms (and higher) in M s and goes to 
zero as Af 0 — » 1 . Also note that the quadrupole and dipole source terms have the same phase. 

Evaluation of Acoustic Pressure 

The previous section describes expressions for the monopole, dipole, and quadrupole source terms 
of LighthilFs equation for shock motion in one dimension. To make statements regarding the relative 
importance of these components to the far field sound, Green’s function is used to solve for the far field 
acoustic pressure. For this analysis, it is assumed that only the first order terms in perturbation Mach 
number M s are necessary to determine the relative importance of the monopole, dipole, and quadrupole 
terms, that the observer is located in the far field, and that the shock is a finite disk in a plug flow 
(fig. 10). The flow outside the plug is ambient and is characterized by the density and the sound 
speed For the analysis presented here, contributions of the interface between the plug flow and 
ambient medium to the far field sound are also neglected. 


jjZii 

dx 1 2 


P~ 



Figure 10. Shock disk in plug flow. 


Pressure From Monopole Source 

To evaluate the far field sound that results from the monopole source, a cylindrical coordinate sys- 
tem is set up with the origin at the shock disk center (fig. 11). The position of the source point along the 
disk is denoted by y and the position to the observer is denoted by x. The magnitude of the vector 
connecting the source point and the observer point R = x - y is denoted by R and is determined by 
geometrical arguments to be 


R = 


i 


r 0 + a 


2 r Q o sin 0 cos(<|> - %) 


(40) 


where r is the radial coordinate, a is the radial coordinate to the source point, <j> is the observer angle 
about the axis of symmetry, and % is the angle to the source point in the plane of the shock disk. 

2 2 

Making the far field approximation o /r 0 —> 0 and employing the binomial expansion, this simpli- 
fies to 


R « r n 


1 - — sin 0 cos(<(> ~X) 
r 0 


(41) 


where r 0 is the distance from the observer to the center of the shock, a is the radial coordinate to the 
source point, 0 is the angle from the jet axis to the observer, <|> is the observer angle about the axis of 
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Figure 1 1. Coordinate system used in evaluation of surface integrals. 


symmetry, and % is the angle to the source point in the plane of the shock disk (fig. 1 1). To solve for the 
far field acoustic pressure from the monopole component of the source term of Lighthill’s equation, the 
source term is multiplied by the free space Green’s function and integrated over the source volume V 
and time x 


P « M = sU, 60p > 


2 


V J x 


m\ + 5 


sin(0)X- /cyj)5(yj -y^)- 




R 


dV dx 


(42) 


where y^ s = 0 is the shock position, and dV = dy^o do d% is the incremental volume. Making the far field 
assumption that 1 !R~ l/r 0 , and integrating overy^. 


— 1 r r ? A/ n M k f R \ 

p mM = jrr L j 60 Pi c oo 2 sin(an-*y^) 8 (x-f- — }dA s dx 


r 0 Ja s J *c f 2 V 

^A/q + 5 J 

where A s is the shock area. Integrating over source time X, 

, , -1 f „ 2 M Q M s k . ( R , 

P M ( X A = 2TT I, 60 Pl c ~ 2 sin I®*- — - k y\s) dA 

0 1 A v c ~ J 

[ M o + 5 J 

-60pici M 0 M s k r « Jf 27t 


4 Jtr n 


M 2 q + 5 


(43) 


r R * r 2n . [ r a sin 0 cos(<b - y)1 1 , , 

sin(cor- r n — ■ — — — \>adody (44) 

J o Jo [ L u Coo-kyi s JJ 


Integrating with the aid of expressions 3.719-2 and 3.715-13 in reference 11, the analytical expres- 
sion for the monopole source pressure is 


P M (x,t) = 


-30p l clkM 0 M s R s 


M 0 + 5 I r, 


(0 sin 


e y i(^ sin6 ^ 


( 


sin 


(Of- 


®V 


V 


(45) 


r 0 
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where /j is the first order Bessel function of the first kind. 


Pressure From Dipole Source 

Again using the Green’s function approach to solve for the pressure due to the dipole component of 
the source term in Lighthill’s equation 


= «/„ I, ~ ~2 *<*-*.> « 

\M Z 0 + 5J 


60p ]Coo MM n 

— — cos(tox -ky x )8'(>- 1 -y l5 ) - 

M 2 0 + 5^ 


(— *.-j 


dV dx 


(46) 


Making the far field approximation MR ~ l/r 0 , integrating the first term over x, and using the fact 
that (ref. 9) 




(47) 


The integral of the first term goes to zero. Making the far field approximation and integrating the 
second term over x yields 


p D (x,t) = - f 

' n « 1 


1 r 60p 


COS 


o J v 


mI + 5 


(48) 


Integrating over V] with equation (47) results in 


P D (*d) 


1 60 P]Caa M s M Q 


r ° ( u 2 ^ * A ’ d)>l 

fid 0 + 5 


L,-|;{ cos K'-£)-^]} 


dA. 


>'l= 0 


1 60p 1 cJW J A/ 0 f . r f R\ if adR\ Jt 
sin co t ky, -k ~ — dA 

h s L V cj C^dyj > 


M 2 0 + 5 


(49) 


where dR/dyi = -cos 0 and dA s = o do d Integrating over the shock area, the expression for the dipole 
source contribution to the far field sound is 


P D (*d) = 


M 0 + 5 | r, 


| & — cos 0 j sin 

CO 

f r \ 

t--* 

c 

^ y 

r R s co sin 0^ 

co sin 0 1 

c 

oo 

_ 

\ °°) 


K 00 ) 


(50) 


r 0 


Note that 


/ cocosG^ 

p d M = p M (x,0(J - — — J 


( 51 ) 
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Pressure From Quadrupole Source 

Using Green’s function, the contribution of the downstream part of the quadrupole source term to 
the far field sound is 


P Q (X> 0 = J^J 10p,c^ 


2 2 AVM M 0-^ 8 


( R 

T -t + — 

V c 


^0 + 5 ) 


R 


cos(cdt — ky j ) dV dx 


(52) 


f — 2 2 2 

where /? = A jr 0 + y 1 4 * a - 2r 0 y 1 cos 0 - 2r Q o sin 0 cos %. Making the far field approximation that 
a / a*q = 0 and y\tr^ ~ 0 and using the binomial expansion R ~ r 0 [l - (yi/r 0 ) cos ® ™ (o/r 0 ) sin 0 cos %]. 
Note that in integrating the quadrupole term, the integral is over the entire source volume since there are 
no delta functions that limit the integral to the surface. The integration over the entire source volume 
provides a challenge for this model because integrating over an infinite volume does not correspond to 
physical jet flows. Therefore, for this analysis, the volume term is limited to a length -4 R s < y\ < 4 R s 
along the axial coordinate. This range is selected to approximately correspond to average shock spacing 
in supersonic nozzles (ref. 12). Therefore, for integration of the downstream quadrupole term, the inte- 
gration limits are 0 < yj < 4 R s . Attempts at obtaining an analytical expression for this integral were 
unsuccessful. Therefore, numerical integration is used to evaluate the quadrupole term. 


To simplify the comparisons and to isolate the sound generated downstream of the sound field, only 
the region downstream of the shock is considered in the volume integration. For the purposes of deter- 
mining the relative importance of the monopole, dipole, and quadrupole terms it is also assumed that the 
source is compact. That is, the wavelength of the sound is large compared with the source size. This use 
of the compact assumption allows an analytical expression for the quadrupole term to be obtained, since 
R in the argument of the cosine term in equation (52) reduces to r 0 . 


Employing the compact source assumption, the analytical expressions for the monopole, dipole, and 
quadrupole terms are found to be 


Pm( x ^ 


-15p ]C ikA s M Q M s 
( M 0 + 5 ) r Q n 


( 


sin 


torn 


(0/ - 


P D (*d) 



to cos 0A 
cjk. J 


(53) 


(54) 


Pq(*,0 


-5o .c^kA M r,M v ? 

(Af 0 - 1 ) sin (kR s ) cos 


( 


M n + 5 r n n 


COr n 


CGf- 


kR 


(55) 


Figure 12 shows the normalized acoustic root-mean-square pressure of the monopole, dipole, and 
quadrupole terms for an observer located 100 shock disk radii (r 0 = 100) away from the mean shock 
position at an angle of 0 = tt/ 4 from the jet axis. For ease of comparison, the root-mean-square pressure 
is normalized by the constant multiplicative factor of the monopole, dipole, and quadrupole pressure 
terms 


l5 Pl clkA s M 0 M s 
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Figure 12. Far field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s analysis of shock noise for three 
frequencies. 


The flow variables are normalized with respect to the sound speed the shock disk radius R s , and 
the upstream density pj. For these comparisons, the wave number k varies with the mean flow Mach 
number as described in equation (20). There are several features in figure 12. First, the quadrupole term 
increases with mean flow Mach number more rapidly with increasing frequency. The dipole terms 
increase with mean flow Mach number. The variation in the monopole terms is an order of magnitude 
smaller than the variation in the dipole terms. Although difficult to see on this scale, the monopole terms 
peak and then decrease with increasing Mach number. The monopole terms dominate for the low fre- 
quencies over the entire Mach number range and for the high frequency case at low Mach number. For 
the case co = 1 , the quadrupole term begins to dominate at about Mq = 2.8. At low frequencies (for which 
the source is highly compact), the contribution of the quadrupole terms is negligible over the Mach 
number range considered in the analysis. The wavelength of the sound is 


X = 


2n Coa j5(\+M 0 ) 


4 


cq a /5 + m; 


o 


(56) 


Thus, low frequencies correspond to highly compact sources. For example for co = 1, the wave- 
length over the range of 1 < Af 0 < 3 is 1 1 .5 < X < 15. The range of wavelengths is not much larger than 
the source size 4/? 5 = 4, and the compact assumption begins to break down. The range of wavelength for 
a source frequency of co = 0.1 is 1 15 < X < 150 and the compact assumption is valid. 

Figure 1 3 shows the effect of observer position on the root-mean-square pressure of the monopole, 
dipole, and quadrupole terms. The results shown are for an observer located at 100 shock disk radii 
away from the mean shock position, at angles of 0 = 0, tc/ 4, and k/2 from the jet axis. The monopole and 
quadrupole terms are unaffected by the change in observer position. The term p M at all observer posi- 
tions shown changes only modestly about a value of approximately 0.016 and pq increases with mean 
flow Mach number from 0 at M 0 = 1 to approximately 0.003 at M 0 = 3.8. The dipole component is 
affected by the observer position. The dipole term is stronger at observer angles closely aligned to the 
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Figure 13. Field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s analysis of shock noise for three 
observer angles for co = 0. 1 . 


jet axis, which is consistent with directivity patterns typical of dipoles. In addition, the dipole term 
reduces to the monopole term at an observer position of 0 = 7t/2. 

Summary of Chapter 1 

Detailed analysis and direct numerical simulation of the sound generated by realistic supersonic jet 
flows are currently impractical. To gain insight into the sound generating mechanisms of such a compli- 
cated flow structure, the elements considered to be essential to the sound-generation process are mod- 
eled. Four elements were identified: Flow past a wavy wall, vortex and vortex interaction, sound and 
shock interaction, and vortex and shock interaction. Because the focus of this work is on noise genera- 
tion in shocked flows, the last two of these elements will be examined numerically in chapters 3 and 4. 

This section presented an analysis of the source terms for Lighthilfs equation for a model of an 
oscillating shock. The results presented for the acoustic pressure are accurate to first order in the pertur- 
bation Mach number and are valid only for an observer in the far field. The figures presented are valid 
for situations in which the source is compact. The analysis shows that the monopole, dipole, and qua- 
drupole terms are all potentially important in shock noise generation. The term that dominates in a par- 
ticular situation depends upon the observer position, the frequency, and the mean flow Mach number. 
For the cases tested here, the monopole term dominates for low frequencies over the entire range of 
Mach numbers studied. Both the monopole and dipole terms are significantly larger than the quadrupole 
term at low frequency. At high frequencies, the quadrupole term dominates at high Mach number. 

Chapter 2 

Numerical Methods and Issues 

The direct simulation of sound generation in shocked flows is challenging because high accuracy is 
required to resolve the acoustic portion of the solution, while dissipation is required to maintain stability 
at the shock. Because shock-generated sound places high demands on the algorithm, care must be exer- 
cised in the selection of the scheme used in the simulations. Two algorithms are used in this report. 
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The first is MacCormack’s scheme, which is second order accurate in time and space. It is used as a 
baseline in scheme comparisons because of its computational efficiency and simplicity. However, 
Gibbs’ oscillations can occur in the vicinity of the moving shock, even with added artificial dissipation. 
MacCormack’s scheme, as it is implemented in this work, is described in the introduction of the section 
entitled “Numerical Error Generated by a Slowly Moving Shock in a Duct.” 

The essentially nonoscillatory scheme is used in most of the calculations presented in this report. It 
was chosen because it allows for high-order accuracy in smooth regions of the flow while minimizing 
oscillations around discontinuities in the solution by adaptive stenciling. Second, third, fourth, and fifth 
order accurate implementations of the ENO scheme were used during the course of this research. The 
scheme is described in the section entitled “Analysis of Exact Solution.” 

The section “Extracting Acoustics From Aerodynamic Flow” examines the calculation of shocked 
flows. There is an additional oscillatory phenomenon which manifests itself when the shock wave is 
moving slowly relative to the mesh. This additional phenomenon has been observed and described and 
has been quantified for high-order schemes (refs. 13 and 14). Much of the content of this section has 
been published previously by the author, but is included here for completeness (ref. 14). 

A high degree of accuracy is required to perform acoustic calculations. High accuracy can be 
obtained by using many cells (or grid points) in a low-order accurate scheme, or by using fewer cells (or 
grid points) in a higher order accurate scheme. A method of determining which approach is most eco- 
nomical is outlined in the section entitled “Economics of Higher Order Schemes.” The result shows that 
the most economical approach depends on the degree of accuracy required. For acoustic calculations 
where the degree of accuracy required is on the order of 10” 7 , the third order ENO scheme proves to be 
most economical. Hence, this is the scheme used in most of the calculations in chapters 4 and 5. 

The section entitled “Extracting Acoustics From Aerodynamic Flow” presents information about 
the extraction of acoustic phenomena from aerodynamic flows. The separation of acoustics from aero- 
dynamics cannot be performed in the near field, and thus the properties of acoustic waves in the far field 
are given. The implications of these properties on setting up and running a calculation for sound gener- 
ation are then described. Some methods of analyzing the acoustics are provided. 


Algorithms 


MacCormack’s Scheme 


The classical MacCormack’s scheme is employed in this research because it is efficient and has 
been used extensively in computational aeroacoustics (refs. 15, 16, and 17). MacCormack’s scheme is 
described in many textbooks, but for completeness the algorithm is briefly outlined here (refs. 18 
and 19). 

Consider the system of hyperbolic equations 


3U 3F 
dt + dx 


0 


(57) 


where U is a vector of J components, x is the streamwise direction coordinate, t is time, and F is a 
J component flux function of U. A simple example of such a system of equations is the Euler equations. 
MacCormack’s scheme approximates the solution to equation (57) by a two step predictor and corrector 
technique. For each component of the U vector, the predictor step is 
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( 58 ) 
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and the corrector step is 



v" + 0 "*' 



(59, 


where the superscript n represents the /ith time level, the subscript j represents the yth cell, and F repre- 
sents the flux function F evaluated at the predictor value U. 

MacCormack’s scheme is second order accurate in space and time. Analysis of the stability of 
MacCormack’s scheme with discussion about the amplification and dispersion factors can be found in 
textbooks such as references 18 and 19. 


MacCormack’s scheme as presented in equations (58) and (59) experiences difficulty maintaining 
stability when applied to shocked flows. Therefore, additional artificial dissipation is required. In this 
algorithm, artificial viscosity of the Jameson type is added to the right hand side of the predictor step 
(ref. 15). The corrected numerical flux is defined by 
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and 


where 
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(62) 


P, + i-2p, + p,-, 

Pm + 2 Pi + Pi- 


Here p is the pressure, u is the flow velocity, and c is the speed of sound. In the calculations presented 
here, the coefficients are defined =1/4 and = 1/256. 

The pressure term in is generally of second order in smooth regions of flow, the term domi- 
nates, and the artificial dissipation is fourth order. Near discontinuities, however, the pressure term 
reduces to zeroth order, the e< 2 ) term dominates, and the artificial dissipation is second order. 

MacCormack’s scheme has numerical damping and dispersion properties since the stencils are used 
throughout the computation. The numerical properties of the scheme are documented in many textbooks 
(ref. 18). 

Because the stencil that is used in the ENO scheme changes spatially as well as temporally, the 
damping and dispersion properties of the ENO scheme are a function of the particular calculation being 
performed. It is not possible to examine these error properties in the traditional manner. Therefore, dis- 
cussion of the properties will focus on specific cases that were studied in this paper. For the reader inter- 
ested in a more quantitative error analysis of the ENO scheme, a recent study of the difficulties in 
obtaining globally high-order-accurate solutions in the simulation of shock-induced sound is recom- 
mended (ref. 21). 
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ENO Scheme 

The class of ENO schemes is chosen because high accuracy is achieved in smooth portions of the 
flow, while spurious oscillations around flow discontinuities such as shocks are bounded (ref. 22). A 
brief discussion of the properties of ENO schemes is given below, but interested readers are referred to 
references 22 and 23 for details. 

An ENO solution operator E h is rth order accurate in the sense of local truncation error 

E h U n = U n+{ + 0(h r+1 ) (64) 

where U(x) is the sufficiently smooth exact solution, and h is the cell spacing. The distinguishing prop- 
erty of ENO schemes is that spurious oscillations near discontinuities in U are bounded. For the one- 
dimensional scalar case, this can be written 

TV(E h U) = TV(U) + 0(h l+q ) (65) 

for some q > 0, where TV{U) is the total variation of U as a function of x , defined as TV(U) = X/ \Ui+\ ~ 
t/,1 for a discrete solution to a scalar conservation law. 

Bounding oscillations near discontinuities is accomplished in ENO schemes through an adaptive 
stencil. Because the stencil used for discretizing the differential equations adapts to the solution, 
schemes based on the ENO property may be thought of as adaptive filters or nonlinear algorithms. 

The ENO schemes minimize numerical oscillations around discontinuities by using data from the 
smoothest part of the flow. At each cell, a searching algorithm determines which portion of the sur- 
rounding flow is smoothest. The stencil spanning this portion of the flow is then used to construct a 
higher order accurate, conservative interpolation to determine the variables at the cell interfaces. 

In this particular finite-volume implementation, the interpolation operator is applied to the cell- 
averaged characteristic variables and the accuracy in space and time is third order. Time integration is 
accomplished by a third order, three stage Runge-Kutta scheme, which is discussed by Shu and Osher 
(ref. 23). The algorithm is applied to the conservation law form of the equations so that shocks are cap- 
tured in the computations and no shock fitting is required to enforce the Rankine-Hugoniot jump rela- 
tions across shocks that appear in the solution. 

Stencil biasing parameters . The ENO schemes can achieve high-order accuracy in smooth regions 
and capture shocks without oscillations by adaptive stenciling. As the ENO schemes were originally 
presented in reference 22, the stencil shifts freely with any detection of a numerical gradient. The direc- 
tion of the stencil shift is determined by the magnitudes of the neighboring divided differences. The 
stencil shifts away from the larger differences. However, a loss of accuracy can occur when this freely 
adaptive algorithm is used (refs. 24 and 25). Shu has suggested that the stencil be biased towards a pre- 
ferred stencil, the one that makes the scheme stable in the sense of linear stability analysis (ref. 26). The 
stencil is allowed to shift only when one neighboring difference is larger than the other by some factor. 
This factor will be referred to as the “bias” parameter. 

This biasing can be accomplished by implementing a factor in the search stencil. To explain how 
biasing is used to affect the stencil, let A* denote the operator that yields the jfcth forward difference on a 
stencil of k + 1 cells with a left index i, which is defined recursively by 


A^w-A* l u k = 2,3,... (66) 
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The algorithm begins by setting j\i) = i. This one-cell stencil results in a piecewise constant recon- 
struction that is spatially first order accurate. To choose k= 1 - 1, the two stencils consid- 

ered as candidates are those obtained by annexing a cell to the left or right of the previously determined 
cell. The stencil that is selected is the one in which the Id h difference is smaller in magnitude 


Ml 

J 


(0 


j k (i) - 1 , if 
j k (i), otherwise 


(67) 


k _ k _ 

where A L u and A R u are the kth differences obtained by annexing the cell to the left or right of the pre- 
viously determined stencil, respectively, (o L ,a /? ) = (l,d) or (a,l) for biasing to the left or right, 
respectively, with a> 1. 

Even when a bias parameter is used, there may be a loss of accuracy when all the numerical gradi- 
ents in a region are small; however, some are orders of magnitude larger than others. Atkins has sug- 
gested another parameter, which serves as a threshold, to force the shift to the preferred stencil 
whenever neighboring differences are small, regardless of their relative magnitudes (ref. 27). This 
parameter will be referred to as the threshold parameter. 

The threshold parameter is implemented as 


I V I j l I i" _L I b _L 1 

if |A l W|<£ and \& R u\ < £ then j ( i ) = j s ( i ) (68) 

k 

where £ is a small parameter and j s (i) identifies the stencil obtained by annexing the fcth cell in the lin- 
early stable direction. 

The results presented in this paper used the definitions cited in references 25 and 26. 


ENO flux computation . The ENO schemes of the Godunov type rely upon the solution of the 
Riemann problem to calculate numerical fluxes. Two methods of computing the fluxes across the cell 
interfaces are evaluated in this report. The first method was developed by Roe and the second was 
developed by Osher and Solomon (refs. 28 and 29). 

Consider a system of hyperbolic equations 


3U 9F 
dt 


= 0 


(69) 


where U is a vector of J components, t is time, x is the streamwise direction coordinate, and F is a J 
component differentiable function of U. 

Roe’s approximate Riemann solver determines the change in flux by finding a mean Jacobian 
matrix A which satisfies 

AF = AAV (70) 

where A represents the difference between any two states in solution space. The matrix A is required to 
have a complete set of right eigenvectors and reduce to the exact Jacobian matrix when the states to the 
left and the right are equal A(U,U) = (3F/3U)U. If Xyj, r^\ and 5 wj are the y‘th eigenvalue, the right 
eigenvector, and the inner product of the jth left eigenvector with AU, respectively, equation (70) can be 
written 


AF = X 


~r (j) i u) Swj 


(71) 
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The flux at the cell interface is written 


/ = < 72 > 
j 

where F L and F R are the values of the fluxes to the left and right of the interface, respectively. Roe gives 
the expressions for X, r, and 8 Wj for the Euler equations in reference 28. 

Osher’s approximate Riemann solver computes fluxes in state space rather than physical space. The 
flux difference between the left and right states is written 

r u * 3F 

iF - i„ t sc M < 73 > 

where the integral is evaluated along an arbitrary path T in state space. 

The flux at the cell interface is given by 

f- + < 74 > 

The evaluation of the integral in equation (74) requires knowledge of the states along each subpath 
T(j) and any sonic states that occur. Osher and Solomon solve for these states explicitly for the Euler 
equations (ref. 29). 

Numerical Error Generated by a Slowly Moving Shock in a Duct 

Introduction 

The numerical treatment of unsteady shocks is challenging. In addition to the usual concerns of sta- 
bility and accuracy, there are conflicting requirements with the calculation of the shock region. While it 
is desirable to resolve the shock crisply by minimizing the smearing effect of artificial dissipation at the 
shock, some artificial dissipation is required at the shock to minimize or eliminate the oscillations that 
occur when attempting to resolve a discontinuity on a finite mesh. There is an additional oscillatory 
phenomenon that can manifest itself in the computation of slowly moving shocks. This additional oscil- 
lation that results from the unsteadiness slowly moving shocks is discussed in references 13, 30, and 31 
and is investigated here to determine the nature of the spurious oscillations and the effect that these spu- 
rious oscillations have on computing sound. 

To simplify the analysis and better isolate difficulties in the numerical calculations, only one- 
dimensional and quasi-one-dimensional flows will be treated here. Thus, vorticity waves will not appear 
and emphasis will be on predicting the acoustic and entropy waves. The model problem to be investi- 
gated is a shock moving at a constant velocity in a one-dimensional flow field. 

Spurious oscillations in unsteady computations of slowly moving shocks have been described 
by Woodward and Colella, who observed the oscillations in computations of high-pressure ratio 
(p^P i > 10 5 ) shocks (ref. 30). These oscillations appear when the speed of the shock relative to the 
mesh is small compared to the maximum flow speed at the shock. Woodward and Colella suggested that 
additional numerical dissipation be added to the scheme at the shock and explained that the reason for 
this spurious numerical behavior is that the shock transition layer alternates between being thick and 
thin as it passes through the mesh. 

Roberts has noted that the oscillation phenomenon is not observed in discontinuous solutions of 
scalar equations (ref. 13). Roberts explains the oscillation in terms of the discrete shock structure. He 
shows that among first order flux-difference splitting schemes, Osher’s approximate Riemann solver 
provides the smallest oscillations because the unsteadiness of the numerical shock structure in state 
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space most closely approximates the true shock structure. Roberts observed the spurious oscillation phe- 
nomenon in calculations for shock pressure ratios as low as 1 .2 in calculations that used first order flux- 
difference splitting schemes. 

Lindquist and Giles observed spurious oscillations in computations that used a Jameson style 
Runge-Kutta scheme with blended second- and fourth-difference artificial dissipation and also with a 
van Leer flux vector splitting algorithm for shock pressure ratios 1.5 < p^Px < 2.1 (ref. 31). They 
described the oscillations in terms of the changing shock shape as the shock traverses the computational 
mesh and found, as did Woodward and Colella, that the spurious oscillations could be reduced by 
smearing the shock over more computational cells. 

Spurious oscillations have been observed during the course of this research in computations that 
used the classical MacCormack’s scheme, a high-order accurate essentially nonoscillatory (ENO) 
scheme, and an implementation of Jameson’s Runge-Kutta scheme, which employs a symmetric total 
variation diminishing (TVD) matrix dissipation. Because the interest here is to compute sound gener- 
ated by shocks and shock and fluid interactions, artificial dissipation is not explicitly increased over the 
solution domain because of the deleterious effect on the generation and the propagation of sound in the 
solution. This lack of increase of artificial dissipation has created some difficulties because the spurious 
oscillations are preserved in the high-order accurate flow computations. 

The ENO schemes of Godunov rely upon the Riemann solver for the calculation of numerical 
fluxes. The effect that two Riemann solvers have on the spurious oscillations will be described in this 
section. In addition, because ENO schemes use an adaptive stencil to reduce spurious oscillations, mod- 
ifications of this adaptive procedure will be examined. Finally, the effects of shock pressure ratio, 
Courant number, grid spacing, and shock speed on the amplitude and the frequency of the spurious 
numerical oscillations will be described. 

Analysis of Exact Solution 

The governing equations for the inviscid, compressible flow in a constant area duct are assumed to 
be the one-dimensional Euler equations 


3U 3F 
dt + 5x 


= 0 


(75) 


where U is the vector of conserved variables [p, pw, pe] T , t is time, x is streamwise direction coordinate, 
and F is the flux vector, [p w, p u + p , ( pe + p)u ] where p is density, u is velocity, e is total energy per 
unit mass, and p is pressure. 


Consider these equations along the duct length 0 < x < L for t > 0 with the initial condition 


U(x,0) = 


U i x < x s 

U 2 x > x^ 


(76) 


where the constant states 1 and 2 represent the flow upstream and downstream of a shock, respectively, 
and Xj is the shock position. 

If these equations are integrated along the duct, then 


v f u <fr + F[U(L,0]-F[U(0,r)] = 0 

dtJ o 


( 77 ) 
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One solution is the trivial case, 1)2 = 112. The nontrivial solution to this equation with the initial con- 
dition as described in equation (75) is a flow field with a shock moving at constant velocity u s , which 
satisfies 


F(U 2 )-F(Uj) = U s ( U 2 -Uj) (78) 

Equation (78) follows directly from equation (77) for constant area states on either side of the 
shock. The solution Uj = U 2 also satisfies equation (78). 

Results 

The unsteady, compressible, inviscid flow in the duct is solved by numerical integration of the one- 
dimensional Euler equations. All variables in the supersonic flow field at the inflow boundary are pre- 
scribed. The static pressure at the downstream (subsonic) boundary is prescribed. Variables are normal- 
ized by the duct length, stagnation pressure, and stagnation sound speed. Control over the shock 
velocity is obtained by making a transformation w = u- u s , where u s is the prescribed shock velocity, so 
that a positive shock velocity moves the shock to the left of the computational domain. 

Computations have been performed over a range of shock pressure ratios and shock speeds, but in 
the interest of brevity, only one typical calculation is shown here. The calculations were performed on a 
512-cell grid at a Courant number of 1. Unless noted otherwise, calculations are performed using Roe’s 
flux solver. Both the bias parameter and threshold parameter are “on,” meaning that the stencil is biased 
towards the preferred one and a threshold limit is set. The values of the biasing parameters that were 
used in these calculations (eqs. (68) and (69)) are a = 2 and e = 1CT 3 . Figure 14 illustrates the spurious 
oscillations observed with the third order ENO scheme for a case in which the pressure ratio across the 
shock is 10.33 and the shock is moving to the left at a speed of 0.05. Figure 14 shows the pressure and 
entropy distributions in the duct after the shock has moved 15 percent of the duct length. Entropy is 
measured by the quantity 5 = p! p Y . Although there were no oscillations at the shock in the initial shock 
position, once the shock begins to move spurious waves develop in the flow solution. The oscillation is 
seen primarily in the entropy wave; the pressure wave is relatively unaffected. 

The spurious error is due to the discrete motion of the shock moving through the mesh. When the 
shock is located at a cell interface, it is extremely thin. As it moves through the cell interior, it smears 
out, weakens in strength, and entropy and pressure waves convect downstream. If the shock-passing fre- 
quency is defined as the frequency associated with the shock passing through a cell 

/shock = £ (79) 


Pressure 



Figure 14. Pressure and entropy as functions of distance along duct. Shock speed = 0.05; shock pressure ratio = 10.33. 
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where Ax is the grid spacing and u s is the shock speed relative to the grid. The wavelengths associated 
with the pressure and entropy waves X p and X s are determined by 


Xp /, 


W 2 + a 2 {U2 + Ci2)Ax 


shock 


x. = 


u 2 A x 


s f, 


shock 


(80) 


(81) 


where w 2 and a 2 are the downstream velocity relative to the shock and sound speed, respectively. When 
adequately resolved on the mesh, the spurious pressure and entropy waves measured in the numerical 
computations compare well with the wavelengths X p and X s described by equations (80) and (81). For 
example, the wavelength X s of the entropy wave in figure 14 is computed as X s ~ 0.550 • 0.00195/0.05 ~ 
0.021. Inspection of figure 14 verifies this result. 

A series of calculations for the same shock strength of p^P\ = 10.333 illustrates the effect of shock 
speed on the behavior of the spurious oscillations. These calculations were performed on a 512-cell grid 
at a Courant number of 1.0. Figure 15 shows the entropy distributions in the tube after the shock has 
moved to the left for a normalized time of 3.0 for shock speeds of 0.02, 0.05, and 0.15. For clarity, the 
entropy distributions are offset by a constant value of 0.1 . The entropy upstream of the shock is 1 .0. The 
wavelengths of the perturbations downstream of the shock are consistent with equation (81). The long 
wavelength disturbances at the lowest shock speeds are only slightly damped downstream of the shock, 
while the short wavelength disturbances at high shock speeds are damped very quickly by the dissipa- 
tion in the scheme. This, of course, explains why these disturbances are not seen when the shock speed 
through the mesh is comparable to the flow speed. 

Figure 16 shows the pressure distributions in the tube after the shock has moved to the left for a nor- 
malized time of 3.0 for shock speeds of 0.02, 0.05, and 0.15. For clarity, the pressure distributions are 
offset by a constant value of 0.05. The pressure upstream of the shock is approximately 0.027. Small 
perturbations are visible in the pressure distribution downstream of the shock and the wavelengths of 
these perturbations are consistent with equation (81). The long wavelength disturbances at the lowest 
shock speeds are only slightly damped downstream of the shock, while the short wavelength distur- 
bances at high shock speeds are damped very quickly by the dissipation in the scheme. 
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Figure 15. Entropy as a function of distance along duct and shock speed. Shock speed = 0.02, 0.05, and 0.15; shock pressure 
ratio = 10.33. 
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Figure 16. Pressure as a function of distance along duct and shock speed. Shock speed = 0.02, 0.05, and 0.15; shock pressure 
ratio = 10.33. 
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Figure 17. Effect of shock velocity on spurious entropy. 

Effect of shock speed. Figure 17 summarizes the effect of the shock velocity relative to the grid on 
the maximum amplitude of the spurious entropy. The ratio of the magnitude of maximum zero-to-peak 
entropy error amplitude to the jump in entropy across the shock is plotted as a function of the ratio of the 
shock velocity (normalized by upstream stagnation sound speed) to upstream Mach number. These cal- 
culations were performed for a shock moving to the right at shock speeds from 0 to 5.0 on a 256-cell 
mesh. The spurious entropy amplitude is machine zero when the shock is stationary relative to the grid. 
The amplitude increases to a maximum when 0.005 < u s /M <0.1 and decreases as the ratio of u s /M 
increases further. 

Figure 1 8 summarizes the effect of the shock velocity relative to the grid on the maximum ampli- 
tude of the spurious pressure. The ratio of the magnitude of maximum zero-to-peak pressure error 
amplitude to the shock strength is plotted as a function of the ratio of the shock velocity (normalized by 
upstream stagnation sound speed) to upstream Mach number. These calculations were performed for a 
shock moving to the right at shock speeds from 0 to 5.0 on a 256-cell mesh. The results show that the 
maximum pressure amplitude is relatively independent of the shock speed and that the maximum error 
in the downstream pressure relative to the shock strength is less than 0.15 percent for all cases studied. 
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Figure 18. Effect of shock velocity on spurious pressure. 



Figure 19. Effect of shock velocity on spurious pressure measured approximately 50 cells downstream of shock. 

Figure 19 shows the zero-to-peak pressure error at a location approximately 50 cells downstream of 
the shock. This figure shows that for faster moving shocks, the amplitude of the pressure error is more 
rapidly damped. This rapid damping is consistent with the observations illustrated in figure 16, which 
relates the wavelength of the spurious pressure with the shock speed. 

Because the amplitude of the error introduced by a slowly moving shock is manifest primarily in 
entropy, further discussion will focus on this flow variable in the section entitled “Effect of Stencil 
Biasing Parameters.” 

Effect of shock strength. Figure 17 also illustrates that the magnitude of the maximum spurious 
entropy amplitude is a function of shock pressure ratio. As the shock strength increases, the magnitude 
of entropy error relative to the static entropy jump decreases. In the weak shock case M = 1 . 1 , the spuri- 
ous oscillations are close to 100 percent of the static jump in entropy over a range of shock speeds. This 
large percentage is because for the weak shock cases, the entropy jump across the shock is very small. 
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(c) Biasing off and threshold off. (d) Biasing off and threshold on. 

Figure 20. Effect of stencil-biasing parameter and threshold parameter on stencil for 256 cells. 


The absolute levels of maximum entropy error are actually orders of magnitude higher for the 
higher Mach number flows. For example at u s IM = 0.04, the absolute levels of maximum entropy error 
are 6s rmx = 0.00033, 0.0039, and 0.097 for Mach numbers M = 1.1, 1.3, and 3.0, respectively. 

Effect of Courant number : A sequence of calculations was performed with Courant numbers 
between 0.1 and 1.0 (the stability limit) to determine the effect of Courant number on the spurious 
entropy. Neither the amplitude nor the wavelength of the spurious oscillations is found to be sensitive to 
the Courant number. This lack of sensitivity is consistent with the results reported by Roberts for flux- 
difference splitting schemes (ref. 13). 

Effect of stencil biasing parameters . In this section the effects of the stencil biasing parameters on 
the spurious entropy are investigated. To illustrate the effects of these parameters on the algorithm, fig- 
ure 20 shows a space and time diagram of the stencil that was used in the ENO algorithm for four 
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Figure 21. Entropy as a function of duct distance for combinations of biasing and threshold parameters for 256 cells. 
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Figure 22. Entropy as a function of duct distance for combinations of biasing and threshold parameters for 1024 cells. 


combinations of stencil biasing parameters. These combinations are (a) threshold on and bias on, 
(b) threshold off and bias on, (c) threshold off and bias off, and (d) threshold on and bias off. For clarity, 
the space and time diagrams are limited to the region near the shock and the total number of cells in the 
computations was reduced to 256. The shock is initially located at x = 0.5 and moves with a velocity of 
0.01 to the left. The upstream Mach number is 1.3. For cases (a) and (d), the threshold is on and a cen- 
tered stencil is used for the majority of the computational space. The centered stencil is the linearly sta- 
ble stencil. Cells on either side of the shock use downwind or upwind stencils as appropriate. For 
cases (b) and (c), the threshold parameter is off and there is a great deal more stencil shifting in the 
smooth regions of the flow. Note that the stencils in these two cases follow entropy and acoustic wave 
paths downstream of the shock. The stencil traces for these wave paths have slopes corresponding to 
acoustic and con vec ting wave speeds. 

Effect of mesh spacing. The effect of mesh spacing is illustrated in figures 21 and 22. The results 
presented in figures 21 and 22 are for a shock moving to the left at a speed of u s = 0.01 with an upstream 
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Mach number of 1.3. The difference in the results occurs because of the difference in discretization. A 
256-cell mesh was used to obtain the results of figure 21 , while a 1024-cell mesh was used to obtain the 
results of figure 22. The entropy distributions in figures 21 and 22 are offset by a constant of 0.0075 and 
the numerical values of entropy are removed from the vertical axis for clarity. Refining the mesh 
reduces the magnitude of entropy oscillation and shows the effect of the biasing parameters more dis- 
tinctly. Each combination of biasing parameters has a unique spurious entropy pattern. Case (a) is 
highly oscillatory with multiple frequencies per wavelength of the oscillation. Case (b) has the smallest 
amplitude of entropy peaks. Case (c) has the fewest peaks per period of spurious oscillation, but the 
amplitude of the largest peak is high. Case (d) has large entropy peaks as well as multiple frequencies 
per oscillation. Although all of the results show significant entropy error, the results obtained by biasing 
the stencil and turning the threshold off (case (b)) provide the lowest amplitude of entropy error. 

Another effect of the mesh spacing is the reduction in the wavelength of the spurious entropy. In 
figure 21, the wavelength of spurious entropy computed by equation (81) is X s ~ 0.34. Refining the 
mesh in figure 22 reduces this wavelength to X s ~ 0.086. The number of points per wavelength of the 
spurious entropy wave is the same for both computations, since u 2 /u s is constant. 

Figure 23 shows the effect of Roe and Osher flux solvers on the spurious entropy. The result is 
shown for case (b) because other cases had similar results. The effect of the Osher solver on the spuri- 
ous entropy was to remove one of the peaks in the entropy wave for each oscillation wavelength that 
was seen in the Roe solver results. 

Calculations of supersonic jet noise will often have shocks moving slowly relative to the grid. 
Although the results presented thus far have shown that spurious entropy exists in calculations with 
slowly moving shocks, it has also been noted that spurious pressure waves are very small in amplitude. 



Figure 23. Entropy as a function of distance along duct length for Roe and Osher flux solvers for biasing on and threshold off. 

Economics of Higher Order Schemes 

The desired accuracy in a computational solution determines whether it is more economical to use a 
second order scheme with many grid points or a high-order scheme with fewer grid points. Figure 24 
illustrates how the cost of implementing these algorithms, as measured by CPU seconds per time step, 
varies as a function of the solution error. 

This figure is constructed in the following manner. The quasi-one-dimensional Euler equations 
described in detail in chapter 4 are solved for the nozzle problem. Grid refinement studies for the 
second, third, and fourth order ENO algorithms are performed on a Cray 2 computer. For each succes- 
sively refined spatial discretization, the error in Mach number is computed for a steady isentropic flow 
in a converging-diverging nozzle. In addition, CPU seconds per time step is measured for the solution 
on each mesh. The relationships of error and CPU time as functions of grid spacing are then combined 
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Figure 24. Computational time per step as a function of error for second, third, and fourth order ENO schemes. 


to construct the figure. Because the temporal and spatial accuracy properties are equivalent for these 
algorithms, the trends presented in figure 24 will be similar whether the accuracy computations are per- 
formed on steady or unsteady flows. The error in figure 24 is measured in the Lj norm, defined by 

N 


where M(x ) is the exact value, M(x) is the numerical approximation, and N is the number of discrete 
values in the numerical solution. 

Figure 24 shows that when an error of the order of 10 -1 is acceptable, the second order scheme is 
more efficient than the fourth order scheme by almost an order of magnitude. If an error on the order of 
only KT 6 is acceptable, then the higher order schemes are an order of magnitude more efficient. Note 
that continually increasing the order of accuracy in a scheme does not necessarily result in a significant 
reduction in CPU time, even at very low acceptable levels of error. Figure 24 shows that for the ENO 
scheme, most of the benefit is realized in moving from second to third order in the range of error norms 
considered. 

It should be noted that the results of figure 24 are particular to the quasi-one-dimensional ENO 
algorithms used in this study. Figure 24 will not apply to the implementation of all schemes, because the 
relative cost of obtaining high-order accuracy for different algorithms will vary. However, the proce- 
dure for determining the relationship of computer cost as a function of acceptable error for different 
algorithms will be the same. 

Extracting Acoustics From Aerodynamic Flow 

Acoustics Defined 

To be able to extract acoustics from an aerodynamic flow, it is important to define the characteris- 
tics of the acoustic waves. For the purposes of this paper, acoustic waves have the following 
characteristics: 

1 . Acoustic waves propagate to the far field at the sum of the local velocity and sound speeds. 

2. Acoustic waves are isentropic. 

3. Acoustic waves are small amplitude. 
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It is often useful to consider geometrical spreading laws to find or categorize acoustic waves. For 
example, when the computation includes a region of uniform (or zero) mean flow velocity, acoustic 
pressure decays proportional to 1/r in three-dimensional flows and \/Jr in two-dimensional flows. 
There is no decay in pressure amplitude in one-dimensional flows. (Note that these decay properties are 
valid when wave steepening is negligible. The spreading characteristics of nonlinear waves are 
discussed in ref. 32.) Geometrical spreading relations are more difficult to interpret in complex flow 
fields. 

The characteristics of acoustic waves can be readily measured in the far field where the amplitudes 
of the acoustic disturbances dominate those of the “hydrodynamic” (nonpropagating) disturbances. 
However, the application of these definitions near the source becomes more difficult. In fact, separating 
the acoustic disturbances from the hydrodynamic disturbances in the acoustic near field is not consid- 
ered to be possible because these quantities do not exist independently (ref. 33). Both hydrodynamic 
and acoustic disturbances are heard by a near field observer, but because the hydrodynamic disturbances 
decay more rapidly with distance from the source than the acoustic disturbances, only the acoustic dis- 
turbances impact the far field. 

Requirements for Acoustic Calculations 

The characteristics of acoustic disturbances determine the requirements of the direct numerical sim- 
ulations used. The numerical simulation must be highly accurate, have a large computational domain, 
and have a long time solution. 

The small amplitude of the acoustic waves relative to the underlying mean flow and the large dis- 
tances the acoustic waves travel require that the numerical simulation be high-order accurate and/or use 
a grid fine enough to sufficiently resolve the acoustic waves. The disparity of length and time scales 
typically of importance in acoustics makes accuracy an issue, since small scales are often important and 
cannot afford to be filtered from the computational solution. In addition, boundary conditions must be 
accurately imposed. 

Because acoustic disturbances are defined as the propagating portion of the unsteady pressure field, 
calculation of acoustics requires a large computational domain. The computational domain must be 
large enough that the pressure field can be measured at least one acoustic wavelength away from the 
source. Acoustic wavelengths are often much larger than the source that generate the sound. 

Long time calculations are necessary to compute at least one period of the far field sound if the 
acoustic signal is periodic. For nonperiodic signals, seven to ten cycles of the lowest frequency are 
required for reasonably accurate spectral estimates. (For a good description of accuracy in spectral esti- 
mates, see ref. 34.) 

Summary of Chapter 2 

The computation of slowly moving shock waves produces spurious, numerical entropy. The spuri- 
ous entropy is a function of the algorithm used in the calculation, and as seen by the modifications made 
to the ENO scheme, even slight changes in a basic algorithm can produce marked changes in the struc- 
ture of the spurious entropy. This phenomenon has been observed by the authors in implementing the 
MacCormack’s scheme, ENO schemes, and a matrix dissipation Runge-Kutta scheme. It has also been 
observed by others who were using flux vector splitting schemes (ref. 31), flux-difference splitting 
schemes with Godunov, Roe, and Osher flux solvers (ref. 13), and PPM methods (ref. 30). Spurious 
entropy normalized by the entropy jump across the shock decreases with increasing shock strength and 
increasing shock velocity, but is insensitive to Courant number. The amplitude of spurious entropy per- 
turbations is relatively unaffected by the flux solver used, but the Osher solver reduces the number of 
peaks in the spurious entropy wave form. 
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Because the amplitude of the spurious entropy wave is a function of the shock speed relative to the 
grid, the obvious method of eliminating the spurious entropy is to move the computational grid with the 
shock during the calculation. Although this grid is not difficult in one-dimensional problems, it is con- 
sidered unreasonable for the multidimensional problems of practical interest, and was not implemented 
during the course of this research. Another approach to reducing the spurious entropy is to increase the 
dissipation of the algorithm, as suggested by Woodward and Colella (ref. 30). This increase in the dissi- 
pation was not implemented, because the added dissipation would affect the acoustic waves as well as 
the entropy. Acoustic pressure waves seem to be less sensitive to numerical errors generated by slowly 
moving shocks. 

In the section “Economics of Higher Order Schemes,” it was shown that the trade-off between the 
added cost of using higher order algorithms depends on the level of accuracy required. The analysis per- 
formed for sound in a converging-diverging nozzle showed that for a numerical error on the order of 
10“*, the third order ENO algorithm is most cost-effective. Therefore, the third order ENO algorithm is 
used for most of the work presented here. 

Chapter 3 

Interaction of Sound With a Planar Shock Wave 

In this section the interaction of a sound wave with a shock is considered. This study is meaningful 
because it models the planar oscillation of a shock wave and shock waves in real aerodynamic flows are 
inherently unsteady. 

In the first three sections, shock capturing formulations of MacCormack’s and ENO schemes are 
applied to the governing equations of fluid dynamics for time dependent, shocked flow through a noz- 
zle. These computations predict the amplification of sound by a shock and compare the solutions pro- 
vided by the different algorithms. Validation of numerical results is always important, but is particularly 
important for these calculations because, although high-order schemes ensure high-order accuracy in 
smooth regions of the solution, they necessarily reduce to first order at the shock. Hence, accuracy is 
lost in the vicinity of the most important region of the solution: the sound source. Since a linear theory 
exists for this problem, numerical results are validated for the case of small amplitude acoustic waves 
incident on the shock. Much of the content in the first three sections was originally published by the 
author (ref. 35), but is included here for completeness. 

The section entitled “Energy Analysis” includes an analysis of the disturbance energy associated 
with a sound wave passing through a shock. This analysis provides insight into the source of distur- 
bance energy generated at the shock. 

Analysis 

Linear Theory 

Within the context of linear theory, only entropy and acoustic waves may exist in a quasi-one- 
dimensional flow (ref. 35). In a fully linear flow, these waves are independent. However, when the flow 
field contains nonlinear features such as shock waves, the nonlinearity acts as a coupling mechanism 
between the linear waves. Thus, the presence of shock waves in a flow field makes it possible for a 
sound wave incident upon a shock to suddenly change its amplitude and generate an entropy wave. 

Entropy is generated across the shock. For a steady flow, this change in entropy is independent of 
time. However, when a sound wave impinges upon the shock, the shock begins to oscillate, the change 
in entropy is no longer constant, and a periodic entropy wave is generated, which convects downstream 
at the local flow velocity. In addition, the impinging sound wave is amplified. 
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Linearized analyses of the interaction of small disturbances with shock waves have been made inde- 
pendently by Blokhintsev (ref. 36), Burgers (ref. 37), Moore (ref. 38), Chang (ref. 39), Kerrebrock 
(ref. 40), and Powell (ref. 41). A numerical study of the interactions of linear plane waves with shocks 
confirms the validity of the linear theory except possibly for waves with incidence angles near or 
beyond the critical angle (ref. 42). 

Landau and Lifshitz report that the ratio of transmitted to incident sound waves determined by lin- 
ear theory is 


(Y - 1 + 2 j 

where A/j is the Mach number upstream of the shock (henceforth called the preshock Mach number), 
A# 2 is the Mach number downstream of the shock, and y is the ratio of specific heats of the fluid 
(ref. 43). Equation (83) and an expression for the ratio of static pressures p^P\ ^ plotted with the ratio 
of specific heats y = 1.4 in figure 25. Note that equation (83) predicts an amplification of the acoustic 
signal as it propagates through the shock wave for all preshock Mach numbers. This amplification is not 
surprising, since the mean flow pressure also increases across a shock. Note, however, that the ratio of 
the perturbation pressures 5p 2 /$Pi is not the same as that for static pressures p^Pv 




2(y - 1 )M^m\ - 1 )- (Af 2 + 1 ) 


(Y 


- 1 )^ + 2 ) 



Figure 25. Ratios of static and perturbation pressures as functions of preshock Mach number. 


Riemann Analysis 

The results of the previous section entitled “Linear Theory” are based on linear theory in which dis- 
turbance amplitudes are assumed to be small. To determine the range for which the linear result is valid, 
the ratio of disturbance pressures across the shock may be computed by considering the iterative solu- 
tion to the Riemann problem. The Riemann analysis is performed by considering the space and time 
diagram illustrated in figure 26. Initially a steady shock wave separates states 1 and 4. At some time A f, 
a disturbance is introduced upstream of the shock which moves the shock and produces acoustic and 
entropy waves downstream. Knowing the initial states 1 and 4, the incident perturbation amplitude, and 
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Figure 26. Quasi-steady Riemann problem analysis for sound and shock interaction. 



Figure 27. Pressure perturbation as function of upstream Mach number. 


utilizing the facts that acoustic waves are isentropic and entropy waves introduce no pressure perturba- 
tion, the ratio of the perturbed states 2 and 3 may be found by numerical iteration. For an incident sound 
wave of pressure perturbation = tp x sin cot, the ratio 8p 2 /8/>| is of interest. The ratio 5p 2 /§Pi is com- 
pared with the linear theory result in figure 27 for several perturbation amplitudes and shows excellent 
agreement for perturbation amplitudes less than e = 1CT 1 . (Results for e < 1CT 2 are visually indistin- 
guishable from the results of Blokhintsev.) Even for perturbation amplitudes of e = 1.0, there is only a 
10-percent difference between solutions at M = 3. 

Combining the expression for 8p 2 /8pj with the Rankine-Hugoniot shock jump relations, similar 
expressions for the fluctuations in density and entropy downstream, as well as for the shock velocity, 
can be calculated. These relations are plotted in figures 28 and 28. Figure 28 shows the relationship 
between the ratios of perturbation pressure, entropy, and density perturbations across the shock. These 
ratios were computed using a perturbation amplitude of E = 10 . This figure shows that the pressure 
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Figure 28. Perturbation ratios as a function of upstream Mach number. 



Figure 29. Shock speed number as a function of upstream Mach number. 


fluctuation is significant over a wide range of Mach numbers. The entropy fluctuations downstream of 
the shock are very small; the quantity 8 5 /A 5 is on the order of e, except near M = 1 where A s — » 0. 
Figure 29 shows the shock Mach number as a function of the upstream Mach number for several pertur- 
bation amplitudes. As expected, the results obtained for small perturbations (e < 1CT 1 ) are indistinguish- 
able from the results of the Riemann analysis, while there is a significant difference in the results 
obtained for the large perturbation (e = 1.0). The shock motion increases with upstream Mach number 
and perturbation amplitude. 

Model Problem 


Governing Equations 

The equations governing the unsteady quasi-one-dimensional flow in a nozzle are 


3U 3F 
dt 


= Q 


(84) 


T 

where U is the vector of conserved variables [p, p u, pe]\ t is time, x is the streamwise direction 
coordinate, Q is a source vector due to the area variation [0, p dA/dx , 0]^, and F is the flux vector 
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9 T 

[p uA, p u + p)A, (pe + p)uA ] , where A is the area of the nozzle, p is density, u is velocity, e is total 
energy per unit mass, and p is pressure. 


Nozzle Shape 

Consider a nozzle which is converging-diverging and designed for a linear Mach number distribu- 
tion when flow is isentropic and fully expanded to produce supersonic velocity in the diverging section. 
The area distribution for such a nozzle is 


Mx) 


l + l—l M (x) 2 
y+ 1 


Y+l 

2 ( 7 - 1 ) 


1 

M(x) 


0 < jc < 1 


(85) 


For the computations performed here, the Mach number at the nozzle inlet is M = 0.8 and varies as 
a function of distance along the nozzle length 

M(x) = (ou + 0.8) 0 < jc < 1 (86) 

where x is the distance along the nozzle, normalized by the nozzle length. 

To simplify the application of boundary conditions at the inflow and outflow, the area distribution 
outside of the nozzle is held constant, that is 

A(x) = A(;c = 0) jc<0 (87) 

A(x) = A(x= 1) jc > 0 (88) 

Two Mach slopes are used in consideration of a range of practical significant preshock Mach num- 
bers to keep the shock close to the center of the nozzle diffuser. The values of Mach slope a used in the 
calculations presented here and the corresponding preshock Mach numbers are presented in table 1. 

A sketch of a nozzle with a Mach slope of 1 is shown in figure 30. 


Table 1. Mach Slopes and Range of Preshock Mach Numbers 


Mach slope 

Pre shock Mach number 

1 

1.4 to 1.8 

2 

2.0 to 2.6 


Governing Equations and Boundary Conditions 

The unsteady, compressible, inviscid flow in a nozzle of varying cross-sectional area is governed by 
the quasi-one-dimensional Euler equations (84). These equations are solved in conservation form so that 
shocks are automatically captured. Total pressure and entropy are prescribed upstream of the shock at 
the nozzle inlet and static pressure is prescribed downstream of the shock at the nozzle exit plane. Flow 
quantities are normalized by the upstream stagnation conditions and the nozzle length L. Finite wave 
conditions developed by Atkins and Casper are employed at the boundaries (ref. 44). These conditions 
are described in appendix C. 

The first step in numerically modeling the interaction of a sound wave with a shock is to determine 
an accurate steady flow solution throughout the nozzle. To obtain sufficiently accurate unsteady results, 
the steady state residual is converged to several orders of magnitude smaller than the smallest perturba- 
tion amplitude to be investigated. Steady flows with residuals only one or two magnitudes smaller than 
the perturbation amplitude may introduce spurious entropy at the inflow boundary in the unsteady 
calculation. 
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Figure 30. Nozzle geometry for Mach slope = 1 . 

To induce unsteady flow through the nozzle, the inflow boundary condition is perturbed sinusoi- 
dally and isentropically. The pressure perturbation is prescribed as 

8p = e/7 sin cor (89) 

where e and co are the normalized amplitude and frequency of the perturbation, and p and t are the 
normalized pressure at the inflow boundary and time. The pressure perturbation introduces an acoustic 
wave at the nozzle inlet which propagates downstream at a velocity equal to the sum of the local veloc- 
ity and sound speed. A range of perturbation amplitudes is selected so that linear as well as nonlinear 
waves can be investigated. The time dependent features of the flow are computed, and the effect of the 
shock on the sound wave is observed. 

Algorithms 

To admit discontinuous solutions, shock capturing formulations of MacCormack’s and ENO algo- 
rithms described in chapter 2 are employed for the computations. Therefore, shock fitting methods that 
explicitly invoke the Rankine-Hugoniot jump relations across a shock are not required. 

Results of Calculations 

Unsteady Calculations 

The results presented in this section are for the case of perturbation amplitude e = 0.01 and fre- 
quency co = 60, which corresponds to an approximate wave number of 5.25. The calculations are per- 
formed on a grid of 128 cells, which corresponds to approximately 24 points per wavelength. The 
nozzle back pressure is prescribed so that a shock appears at x = 0.6 in the steady solution. For a nozzle 
with a Mach slope a = 1 , this condition corresponds to a preshock Mach number of 1 .4. 

Figure 31 shows the perturbation pressure, density, and velocity in the nozzle at a normalized time 
of n and the results provided by MacCormack’s and the second, third, and fourth order ENO schemes. 
Perturbation quantities are determined by subtracting the mean flow quantities from their time depen- 
dent counterparts 

Sp(xj) = p{xj) -p(x 9 0) 

5p(x,/) = p(x,t) - p(*,0) 

8 u(x,t) = u(x,t) - w(jc,0) (90) 
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MacCormack's 

Second order ENO 

Third order ENO 

• Fourth order ENO 


Figure 31. Pressure, density, and velocity perturbations along nozzle length. 


where p is pressure, p is density, u is velocity, x is streamwise directional coordinate, and t is time. 

The MacCormack’s and the third and fourth order ENO schemes do an excellent job of predicting 
the perturbation pressure upstream of the shock where the flow is smooth. The second order ENO 
scheme, which behaves similar to a TVD scheme because it has only second order interpolation, has a 
slight leading phase error. At the shock, the ENO schemes do better at capturing the shock in fewer 
cells. The flow solutions differ more significantly downstream of the shock. The phase shift error is 
amplified and the wave amplitude is dissipated in the second order ENO results. The MacCormack 
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Figure 32. Pressure perturbation along nozzle, e = 10 -5 ; Afj = 1.58. 



Figure 33. Pressure perturbation along nozzle. E = 10 5 ; = 2.36. 


solution is slightly damped relative to the third and fourth order ENO solutions and has spurious oscilla- 
tions on portions of the perturbations close to the shock. The difference in the solution between the third 
and fourth order ENO schemes is not graphically perceptible on this scale. Note that with mesh refine- 
ment all schemes (including MacCormack’s and the second order ENO scheme) converge to the solu- 
tions produced by the higher order ENO schemes. 

These results show that all the algorithms applied to this problem compute the perturbation quanti- 
ties in the nozzle well. The third and fourth order ENO schemes provide the best solutions because they 
have no oscillations at the shock and less damping downstream. The cost of running high-order ENO is 
higher than second order ENO or MacCormack’s scheme, however. The economics of high-order 
schemes is discussed in the section entitled “Economics of Higher Order Schemes.” 

Effect of Mach Number 

Sound amplification by a shock wave is highly dependent upon the preshock Mach number. Fig- 
ures 32 and 33 show a sequence of pressure perturbations moving through the nozzle for preshock Mach 
numbers of 1 .58 and 2.36. A small amplitude perturbation of e = 10 -5 is introduced at the inflow bound- 
ary. The perturbation maintains its sinusoidal shape throughout the nozzle. The amplitude decreases as 
the flow expands in the nozzle upstream of the shock. At the shock, it is amplified as the flow is 
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compressed and continues to increase gradually as the flow is compressed further. The amplification of 
the sound wave is much more significant at higher Mach number. 

Comparison of Numerical Results With Linear Theory 

To quantify the effect of Mach number on the perturbation amplification and validate the computa- 
tions, the numerically determined pressure perturbation amplitude ratio is compared in figure 34 with 
the linear theory over a range of Mach numbers for perturbation amplitudes of e = 10 -5 and e = 10 -2 . 
The first amplitude is small enough to be well within the validity of the linear theory. The second ampli- 
tude approaches the limit of linear theory validity, particularly for lower Mach numbers where the shock 
is weak and may undergo large excursions. The normalized acoustic perturbation wave number at the 
inflow boundary is set to 4. Results for Mach numbers between 1.4 and 2.6 are presented here. The 
computations have from 22 to 27 cells per wavelength. Table 2 lists the cells per wavelength for the 
cases presented in this paper. 

Figure 34 compares the numerical results with linear theory. The third order ENO scheme and lin- 
ear theory match extremely well for the very small amplitude perturbation of e = 10~ 5 . The differences 
between numerical solutions and linear theory become more pronounced at the higher perturbation 
amplitude of e = 10 -2 . Some discrepancy between the linear theory and numerical results is not surpris- 
ing since nonlinearities may become important at this perturbation amplitude. Because the shock in 
the MacCormack solution is spread over several cells, determination of the ratio 5/? 2 /§/?] is more diffi- 
cult, particularly when the shock is located close to the end of the nozzle. This difficulty is why no 
MacCormack result is shown for M = 2.6. 


Table 2. Preshock Mach Number and Minimum Cells 
per Wavelength for Calculations 


Preshock Mach 
number 

Minimum cells per 
wavelength 

1.4 

26.92 

1.58 

26.18 

1.7 

25.75 

2.0 

23.52 

2.36 

22.50 

2.6 

22.13 



Figure 34. Perturbation ratio as a function of preshock Mach number. 


45 






Energy Analysis 

The previous sections have shown that according to linear theory and numerical computation, sound 
pressure is amplified as a sound wave passes through a shock. This section will address the issue of 
whether acoustic energy is generated in the interaction process. 

To address this issue, consider Myers’ energy corollary (ref. 45). Myers’ energy corollary is an 
exact equation governing the transport of energy associated with an arbitrary flow field. This corollary 
describes the acoustic energy in general situations where the linear description of energy is inadequate. 
Clearly, with the highly nonlinear flow field associated with sound and shock interaction, such an 
approach is warranted here. 

For a one-dimensional flow in which viscous effects are negligible, Myers’ corollary reduces to 


dE dW 
dt 


= -D 


(91) 


where E is energy density, t is time, W is disturbance energy flux, and x is streamwise direction coordi- 
nate and the disturbance energy density divided by area is defined 


E \ f 2 2\ 

A = 2 P V - “oJ + p o“o ■ (“0 “ - *o) - (/> - Po) + P(A - *o) ( 9 2 ) 

where A is the cross-sectional area of the duct, T is the absolute temperature, h is the enthalpy, 5 is 
entropy, p is density, p is pressure, and u is velocity. The subscript 0 represents the undisturbed state. 
The disturbance energy flux divided by area is 

W 

— = (P«-p 0 “oH*-^o“ 7 \)( 5 -*o)] + Po w o( 7, ~ 7 V(‘ y ~*o) ( 93 ) 

and the disturbance energy source divided by area is 

J = -(PM-P 0 H 0 )(5 — J 0 ) • Vr 0 -(i-J 0 )p o ii 0 ' V(7’-7 0 ) (94) 

For the one-dimensional sound and shock interaction, only acoustic and entropy modes are present 
(vorticity requires three dimensions). Thus, the energy density can be divided into two types of energy, 
acoustic E a and entropy E s , as 



dh 

p 0 u Q ( u Q -u)-(p-p 0 ) + p^ 


(p-p o) 


(95) 


and 
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( 97 ) 
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and 




( 98 ) 


where y is the ratio of specific heats for an ideal gas. 

It is instructive to examine the energy components and energy source terms in space and time to see 
how these quantities change during the sound and shock interaction process. These quantities from 
equations (95), (96), and (94) are shown in figures 35 to 37 for the Mach number upstream of the shock 
M\ = 3, the acoustic disturbance amplitude e = 0. 1, and 512 cells, which are distributed along the duct 
length. For these calculations, the shock is initially at r = 0.4. The sound disturbance enters the duct at a 
time of t = 0 at the duct inlet and propagates downstream. At a time of f « 0.13 the sound wave hits the 
shock. After the interaction, figure 35 shows that acoustic energy is present downstream. The amplitude 
of the sound energy downstream of the shock is higher than the energy upstream, which indicates that 
acoustic energy is generated by the sound and shock interaction process. Note that the slope of the path 
of the sound wave in space and time increases after interaction with the shock. The inverse of the slopes 
before and after the shock corresponds to the quantities Uy + C\ and m 2 + c 2 , respectively, where u is 
velocity, c is speed of sound, and the subscripts 1 and 2 refer to the states upstream and downstream of 
the shock, respectively. 

Figure 36 shows the component of entropy energy in space and time. Because the sound wave is by 
definition isentropic, there is no disturbance entropy upstream of the shock. After the sound wave hits 
the shock at t « 0.13, however, entropy energy appears downstream. Thus, entropy energy is generated 
during sound and shock interaction. The inverse slope of the path of the disturbance entropy corre- 
sponds to the downstream convection velocity, as expected. Figure 37 shows the variation in the source 
term of equation (94) and provides insight into the generation of disturbance acoustic and entropy 
energy downstream of the shock. The source term is zero except along the shock wave for t > 0.13. 
Thus, the shock wave is the source of the disturbance energies downstream. Equation (94) indicates that 
the source of disturbance energy is the transfer of energy from the mean flow. This transfer can occur in 
the presence of a temperature gradient with fluctuations in entropy and momentum and with 
fluctuations in entropy and temperature. Note that the source term is zero when the fluctuation in 
entropy is zero. This correlation implies that the disturbance energy source is related to the shock 
motion. 



Acoustic 

energy 


0.068 

0.058 

0.049 

0.039 

0.029 

0.019 

0.009 

-O.(KX) 

- 0.010 

- 0.020 

-0.030 


x 


Figure 35. Disturbance acoustic energy as a function of space time for 512 cells distributed along duct length. Preshock Mach 
number = 3; disturbance acoustic amplitude = 0.1. 
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Figure 36. Disturbance entropy energy as a function of space time for 512 cells distributed along duct length. Preshock Mach 
number = 3; disturbance acoustic amplitude = 0.1. 
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Figure 37. Disturbance energy source as a function of space time for 512 cells distributed along duct length. Preshock Mach 
number = 3; disturbance acoustic amplitude =0.1. 


Summary of Chapter 3 

This section describes a numerical investigation of sound amplification by a shock wave. 
MacCormack’s and the second, third, and fourth order ENO schemes are employed to compute the time 
dependent shocked flow field in a converging-diverging nozzle. The flow is disturbed by introducing a 
sound wave at the inflow boundary. All schemes are shown to predict the perturbation amplitude and 
phase speed in the nozzle well, especially in the supersonic, smooth portion of the flow. The high-order 
ENO schemes provide the best overall results because the flow around the shock does not contain spuri- 
ous oscillations and the flow downstream of the shock shows little dissipation. 

The numerical results are compared with linear theory. Linear theory validates the numerical solu- 
tion for small perturbation amplitudes. Numerical results for the larger perturbation amplitude also com- 
pare well with the linear theory, which indicates that the linear theory has a wide range of applicability 
in predicting sound amplification by a shock. 

Analysis of the equation governing the perturbation energy shows that disturbance energy is gener- 
ated at the shock. The source term for this energy goes to zero as the entropy fluctuation goes to zero. 
This correlation indicates the importance of shock motion on the generation of disturbance energy. 
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Chapter 4 


Interaction of a Vortex Ring With a Shock Wave 

In this section the interaction of a vortex ring with a normal shock wave is considered. The study of 
this interaction models a fundamental mechanism of sound generation in supersonic jets: the interaction 
of turbulence with shock waves. 

Early research on shock and vortex interaction focused primarily on experimental studies (refs. 46, 
47, and 48) and the development of predictive linear theories (refs. 49, 50, and 38), which were com- 
pared with experimental results. Ribner suggested that the study of the interaction of a vorticity wave 
with a shock would provide useful information about the sound generated by turbulence in supersonic 
jet flows (ref. 49). Ribner studied the problem analytically and developed a linear theory that described 
sound generated by shock-vortex interaction (refs. 49 and 50). Moore studied the interactions of a vari- 
ety of plane wave disturbances with shocks in an unsteady reference frame (ref. 38). 

Pao and Salas studied the interaction of a shock wave and a vortex numerically (ref. 51). Their study 
investigated the interaction of a columnar vortex with a normal shock wave by solving the Euler equa- 
tions and using MacCormack’s scheme with a shock-fitting numerical technique. Salas (ref. 52), 
Hussanini et al. (ref. 53), and Kopriva et al. (ref. 54) applied spectral methods with shock fitting meth- 
ods to the shock-vortex problem. Spectral methods provided increased accuracy of the solution, but 
were limited to weak shock and vortex interaction cases. Meadows, Kumar, and Hussaini (ref. 55) stud- 
ied the interaction of a columnar vortex with a shock wave by using a shock capturing scheme. Shock 
capturing proved to be beneficial because strong shock and vortex interaction cases that result in the for- 
mation of secondary shocks could be studied readily. The authors noted that to provide a quantitatively 
accurate representation of the acoustic wave, improved downstream boundary conditions and higher 
order numerical schemes were required. Casper then investigated the shock and vortex interaction prob- 
lem with a high-order ENO scheme and found that higher accuracy greatly improved the resolution of 
the acoustic wave downstream of the shock (ref. 56). 

Reference 56 studied columnar vortex and shock interactions. The interaction of a vortex ring with 
a shock closely models the interaction of turbulence within the shear layer with shock waves present in 
the plumes of imperfectly expanded axisymmetric supersonic jets. 

Model 

The calculations presented in this section model the interaction of a vortical structure within the 
supersonic portion of the shear layer with a shock wave. The core size is small relative to the size of the 
ring in these calculations, primarily because of numerical considerations. Because the shear layer in a 
jet spreads with increasing distance from the nozzle exit plane, a vortex ring with a small core models 
the interaction of disturbances in the shock cell closest to the nozzle lip. When results are presented in 
dimensional units, the variables are dimensionalized with ambient stagnation sound speed and atmo- 
spheric pressure; thus, the calculations closely model disturbances that interact with the first shock cell 
of an overexpanded nozzle (fig. 37). 

The strength of the vortex is chosen to correspond to the observed strengths of turbulent fluctuations 
in the jet shear layer. Experimental observations indicate that an appropriate level of velocity distur- 
bance is approximately 3 percent of the mean flow velocity. The vortex strength chosen for the majority 
of the calculations presented here (A/, = 1 .5) has a velocity perturbation (the ratio of the vortex core 
velocity to the upstream mean flow velocity) of 6.8 percent. 

Although imperfectly expanded jet flows typically contain systems of oblique shocks, this 
model uses the interaction of a ring vortex with a normal shock wave. This model is a reasonable 
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Figure 37. Model of vortex ring and shock interaction. 

approximation because the oblique shocks tend to curve upstream with increasing distance from the jet 
centerline and are approximately normal upon termination in the shear layer (ref. 58). 

Vortex rotation is taken to be counterclockwise for the majority of the calculations because the 
velocity decreases with increasing distance from the jet axis, which results in positive vorticity. Some 
results are shown for a clockwise rotating vortex in the section entitled “Typical Interaction of Clock- 
wise Vortex/' 

Geometry 

The geometric model of the interaction of a ring vortex with a shock wave that was used in the com- 
putations is illustrated in figure 38. A vortex ring is introduced upstream of the shock wave. The vortex 
ring is characterized by its strength T, its core radius r c , and the distance from the axis of symmetry to 
the vortex filament r 0 . The shock is characterized by the Mach number of the flow upstream of the 
shock M j . 

Governing Equations 

The equations that govern the interaction of a vortex ring with a shock wave, neglecting viscous 
effects, are the Euler equations of gas dynamics in axisymmetric coordinates: 


where 


F = 

G = 

H = [0,0,p,0] r (103) 


drQ drF drG u 

~df + + = H 


(99) 


Q = [p,p«,pv,p<?] 

2 T 
Ip«,pw + p,puv,(pe + p)u] 

2 T 

[pv,pv«,pv + p,(pe + p)v] 


( 100 ) 
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Figure 38. Vortex ring and shock interaction. 


p is density, u and v are the velocity components in the axial and radial directions, e is total energy per 
unit mass, p is pressure, x is the axial position, and r is the radial distance from the axis of symmetry. 

Boundary Conditions 

Conditions are prescribed at the inflow and outflow computational boundaries to establish a shock 
with the specified strength at x = 7. The boundary condition along the axis of symmetry is prescribed by 
requiring that the fluxes across the axis are zero. The method of prescribing the boundary conditions is 
described in reference 44. 

Solution Procedure 

The calculation is performed in two steps. First, a steady, shocked flow is established with the flow 
parallel to the axis of symmetry. Next, a vortex ring is introduced to the flow field 7 core radii upstream 
of the shock at a time of T = 0. All flow variables are normalized with respect to the static pressure and 
density upstream of the shock and the size of the vortex core radius. The vortex ring then convects 
downstream and passes through the shock wave. After this interaction, significant changes occur down- 
stream of the shock wave. Sound, vorticity, and entropy are generated as the vortex interacts with the 
shock wave. The primary focus of this section is on the acoustic waves generated by the interaction. 


Vortex Model 

A model of the ring vortex introduces appropriate perturbations in velocity and pressure to the flow 
field as an initial condition. The vortex is introduced at the initial time t = 0 and at subsequent times the 
Euler equations determine the entire flow field — including the vortex. At all subsequent time steps the 
Euler equations are satisfied to a level of error that corresponds to the truncation error of the ENO algo- 
rithm, which is third order accurate in space and time. 
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The initial condition prescribed for the vortex is the classical toroidal vortex of Lamb, which is aug- 
mented with a solid body core (ref. 59). The derivation of the velocity and pressure equations is pro- 
vided in appendix B. 

Limits of the Vortex Ring Model 

The cross section of the vortex ring core is approximately circular when the vortex filament is 170 
core radii away from the axis of symmetry (ref. 60). When the vortex core radius increases relative to 
the ring radius, the degree of circularity decreases, until the vortex becomes Hill’s vortex (ref. 60). The 
streamlines of this vortex become oblong and flatten close to the axis of symmetry. To closely approxi- 
mate the vortex ring with a circular core, the vortex filament is located 125 core radii away from the 
axis of symmetry unless otherwise noted. The small size of the vortex core radius relative to the ring 
radius reasonably approximates the size of a turbulent structure within the supersonic portion of the 
shear layer near the first shock cell. 


Vortex Parameter Modeling 

The parameters of the vortex ring which model the physical jet turbulence and shock interaction 
will be described here. The first parameter of interest is the vortex strength. Circulation is typically con- 
sidered to be a measure of the vortex strength. For the vortex ring, the circulation is related to the geo- 
metrical parameters and induced velocity V according to (ref. 60) 


T = 


4V7tr 0 


In — 

v r < J 


1 

4 


(104) 


The next parameter of interest is the ratio of the core radius to the ring radius r c /r 0 . The core radius 
models the size of the disturbance in the supersonic portion of the mixing layer when it interacts with 
the shock wave. The ring radius models the radial distance between the jet axis of symmetry and the 
inner boundary of the turbulent mixing layer. This ratio is a function of distance downstream from the 
jet exit plane. 


Vortex Preservation Study 

An important feature of this calculation is the accurate representation of the source of the sound 
generation: the interaction of a vortex with the shock. To ensure that the vortex is adequately resolved in 
the calculation, the time history of the minimum pressure within the vortex is tracked. A series of com- 
putations is performed where the vortex ring convects over grids of various resolutions. It was found 
that while convecting on a uniform mesh with 10 cells per core diameter, the minimum pressure varies 
by a maximum of 0.035 percent over the time it takes the vortex to travel 7 core radii, which is the initial 
distance of the vortex from the shock. This level of numerical error is considered to be acceptable; 
therefore, this grid resolution was held fixed for all the computations. 


Computational Grid 

As mentioned in the previous section, a grid with 10 cells per vortex core diameter adequately 
resolved the vortex in the computation. This grid resolution is used in the fine uniform grid in the region 
-6 < x < 85 and 72 < r < 178, where Ax = Ar = 0.2. However, it is computationally prohibitive to use a 
uniform mesh with such fine resolution throughout the entire computational domain. The computational 
domain is required to be large enough so that at least one wavelength of the acoustic disturbance is con- 
tained in the calculation. This size of domain allows accurate measurements to be made in the acoustic 
far field so that computations to determine the acoustic energy level can be made. In addition, the com- 
putational domain must also be large enough that once the acoustic wave has passed a far field observer 
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point, the wave has time to go through one period of its oscillation without being contaminated by errors 
introduced by the boundaries. Thus, in the calculations shown here, the computational grid is uniform 
over the range of the acoustic source interaction and sound wave propagation through at least one 
acoustic wavelength. The grid is then stretched to the boundaries by using a hyperbolic sine function. 
The grid contains 481 cells in the axial direction and 556 cells in the radial direction. Every fifteenth 
cell of the grid is shown in figure 39. 

Note that this grid resolves the generated acoustic waves. Preliminary test calculations performed 
over Mach number range 1.1 < M < 1.7 cover an acoustic peak-to- valley wavelength range of 2.5 < X 
< 3.0. Thus, the minimum resolution of the acoustic wave is 12.5 points per peak-to- valley wavelength. 
The actual number of points per wavelength is at least twice this number. Thus, the acoustic wave is 
well resolved. 

Also note that there is another length scale of possible interest in this problem, which is not well 
resolved on this grid. During the interaction of the vortex with the shock, the shock wave begins to 
move. For typically small vortex strengths, the shock displacement is small. In fact, for many of the cal- 
culations performed during the course of this research, the shock moves through the distance of only 
one cell. Because further grid refinement proved to be prohibitively expensive and one-dimensional cal- 
culations of sound and shock interaction provided results which matched linear theory well without 
resolving the shock motion, the grid was not refined further. 

Typical Interaction of Counterclockwise Vortex 

A typical interaction of a counterclockwise rotating vortex with a shock wave is described in this 
section. In the calculations for this case, the Mach number upstream of the shock is M = 1 .5 and the vor- 
tex strength is T = 0.75. 

In the contour plots shown in this chapter, contours are shown for a part of the computational 
domain near the interaction point of the vortex core and the shock wave. Because of the axisymmetry, 
only a single cross section of the solution will be shown in the results which follow. All the results pre- 
sented are within the uniform grid region of the computational grid so that effects of grid stretching are 
not present in the results. The contours are computed on a grid which is five times as coarse as the com- 
putational grid. Note that the contour range is kept constant for each plot so that relative values of the 
perturbation quantities at different time levels can be compared. In these calculations, a unit of time is 
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defined as 7 = r L Ju x where w, is the upstream flow speed; thus 7 represents the time it takes the vortex 
core to move 1 radius in the flow upstream of the shock. 

Contour plots of the flow variables are shown for three selected times: 7 = 0, 7 = 8, and 7 = 50. At 
time 7=0, the vortex is upstream of the shock. At 7 = 8 the most upstream edge of the vortex core is 
approximately aligned with the shock; at 7 = 50, the vortex is approximately 30 core radii downstream 
of the shock. 

Pressure 

Figure 40 shows the change in pressure from the mean state. The only perturbation in pressure at 
7 = 0 is the decrease in pressure at the vortex core. As the vortex begins to interact with the shock, 


Figure 40. Contours of pressure perturbation downstream of shock at T ~ 0. T= S. and T = 50. 
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Figure 41 . Peak pressure perturbation as a function of time. 


additional pressure perturbations are generated downstream of the shock. At T = 8 high amplitude pres- 
sure disturbances are seen just downstream of the shock. As time passes, these pressure disturbances 
travel downstream more rapidly than the vortex. Figure 41 shows a plot of the position of the peak pres- 
sure as a function of time. The plot shows a linearly increasing change in the position of the pressure 
disturbance as a function of time, which means that the pressure disturbance is traveling at a constant 
speed relative to the vortex. The slope of this curve is found to be AR/AT = 0.77 core radii per period, 
which is the sound speed downstream of the shock. (The velocity 0.77 core radii per period is equivalent 
to a velocity normalized by the sound speed upstream of the shock of 1.15. The ratio of the sound 
speeds across the shock is 1.15 for an upstream Mach number of 1.5.) Thus, the pressure disturbances 
are traveling at the sound velocity relative to the mean flow, satisfying one of the defining features of 
acoustic waves. 

The structure of the acoustic wave is readily apparent at T = 50 in figure 40. As predicted in linear 
theory, the sound wave is quadrupole in nature. The acoustic wave front is composed of alternate com- 
pression and rarefaction fronts. 

Pressure perturbations along radii extending from the vortex center at 10° increments are shown in 
figure 42. This figure quantifies the change in amplitude as a function of the angle from the horizontal 
and shows that the peak-to-peak pressure amplitude is maximum at 0 = 50° and 0 = -55°. From this fig- 
ure, the peak-to-valley measure of the wavelength of the pressure disturbances is seen to be X ~ 2.8. 

Figure 43 shows a carpet plot of the pressure levels downstream of the shock at T - 50. The image 
is processed with the Fast computer software (ref. 61) and is shown at an angle to clarify the detailed 
features of the flow field. Note the resolution of the cylindrical acoustic wave and the complex system 
of pressure waves between the shock and the cylindrical acoustic wave front. 

Figure 44 shows the product of the square root of the distance the acoustic wave has traveled (rela- 
tive to the vortex core) and the peak pressure amplitude along a line extending from the vortex core at 
60° from the horizontal passing through the vortex filament as a function of distance from the vortex fil- 
ament. The product varies significantly at small distances from the core, which is synonymous with 
early time and proximity to the source, but flattens to an almost constant level at larger distances. This 
behavior of the product shows that in the far field, the sound pressure decays as 1 1 Jr, which is charac- 
teristic of geometrical spreading in two dimensions (cylindrical spreading). Although these calculations 
are performed in an axisymmetric coordinate system, which should allow for three-dimensional geo- 
metrical spreading of the sound, the presence of the shock wave prohibits the spreading upstream. In 
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Figure 42. Pressure perturbations along radii extending from vortex core at x - 30. 



Figure 43. Pressure perturbations downstream of shock at T= 50. 

realistic flows, where viscosity plays a role in the fluid dynamics, shock waves are finite in extent. In 
these model flows, cylindrical spreading would occur close to the shock wave, while spherical spread- 
ing (with p decaying as 1/r) would occur in the far field. Hardin has shown that in cylindrical coordi- 
nates, the sound decay rate is a function of the source size to the distance from the source to the observer 
(ref. 62). In the section “Shock Dynamics,” it will be shown that the interaction of the vortex ring and 
shock wave produces a disturbance that travels along the shock, thereby increasing the size of the 
potential noise source with increasing time. 

A cusp-shaped pressure wave in the pressure field, which results from the ring vortex and shock 
interaction, is apparent in the flow field downstream of the shock and below the vortex filament. The 
ring vortex and shock interaction produces not only the acoustic quadrupole, but also a cusp-shaped 
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Figure 44. Decay rate of acoustic pressure. 


pressure wave that connects the portion of the acoustic wave immediately below the vortex filament and 
the shock wave. This feature was not found in columnar vortex and shock interaction studies and 
appears to be sensitive to the introduction of the initial condition (refs. 46, 49, 50, 51, 55, and 56). 
Because of computational cost, the vortex is introduced seven core radii upstream of the shock for most 
of the calculations presented in this paper. However, to test the effect of the initial distance between the 
vortex and the shock wave, Ax j on the flow field downstream of the shock, a single calculation is pre- 
sented where the vortex is initially 23 core radii upstream of the shock. Figure 45 shows the solution 
downstream of the shock at times T = 50 and T - 66 for initial shock locations of x = 7 and x = 23, 
respectively. In both of these cases the vortex has traveled approximately 23 core radii after its filament 
has passed through the shock wave. The cylindrical portion of the pressure field is similar for both ini- 
tial conditions; however, there are differences in the downstream flow field. Most notably, the cusp- 
shaped pressure structure is much larger for the longer time solution (Axj = 23), which indicates that it 
originates at the beginning of the computation. Further study of this feature is required to quantify the 
effect of axisymmetry and the initial condition on its structure, and to determine whether it is numerical 
or physical in nature. 

Density 

Figure 46 shows the change in density of the flow field. At time T = 0 the only change in density is 
associated with the vortex. As the vortex core begins to interact with the shock, density waves appear 
downstream of the shock. The interesting character of the density waves becomes clear in the final fig- 
ure of the time sequence. Two types of density perturbations are evident in the figure. Density distur- 
bances associated with the acoustic wave propagate in a nearly circular pattern downstream of the 
shock. This is to be expected, because acoustic waves are by definition isentropic and there is a clear 
relationship between the density and pressure (s = constant = p! p^). In addition, there are also convec- 
tive density disturbances. These density disturbances look like spokes reaching out from the vortex core 
and terminating at the shock. These disturbances are associated with the entropy waves that are associ- 
ated with changes in the shock strength. As the vortex interacts with the shock, the shock wave begins to 
move. As the shock moves, the change in entropy across the shock is no longer constant and an entropy 
wave is generated that convects downstream at the local flow velocity. For the counterclockwise rotat- 
ing vortex, the portion of the shock above the vortex filament initially moves upstream and the portion 
of the shock below the filament moves downstream. As time increases, disturbances move away from 
the interaction point and their motion along the shock wave creates entropy perturbations that convect 
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(a) Vortex initially 7 core radii from the shock at T = 50. 
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(h) Vortex initially 23 core radii from the shock at T = 66. 

Figure 45. Pressure contours downstream. 

downstream at the local flow velocity. These disturbances show a strong resemblance to the features 
observed by Naumann and Hermanns and described as contact surfaces (ref. 46). 

Figure 47 shows a digitized version of a Mach-Zehnder interferogram from reference 46 that illus- 
trates the sound wave and contact surface observed in the experiment. In the experiment, the interaction 
is produced as follows. A sharp edged profile is placed in a shock tube. A diaphragm (located to the left 
of the airfoil) is broken and a weak disturbance travels through the tube and over the trailing edge, 
which produces a starting vortex. Another diaphragm (located at the right of the airfoil) is broken and a 
shock wave travels towards the vortex. The results of the interaction are shown in figure 47. Because the 
interaction observed in this experiment is strong, the motion of the shock wave is pronounced and the 
acoustic wave front is asymmetrical. The contact surfaces form a curved funnel-shaped structure 
between the shock wave and the vortex. 

Figure 48 shows a carpet plot of the density at T = 50. The sound wave and complex nature of the 
contact surfaces are clearly visible. Experiments provide evidence for the physical nature of the contact 
surfaces observed in these computations. However, earlier analysis in chapter 3 demonstrated that the 
computation of slowly moving shock waves produces error that manifests itself primarily in entropy. 
Therefore, caution must be exercised in the interpretation of the strength of these disturbances. Further 
analysis is required before the nature of these contact surfaces is validated. 
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Figure 46. Contours of density perturbation downstream of shock at T - 0, T - 8, and T = 50. 


Vorticity. Figure 49 shows the change in vorticity. At the initial time, the only vorticity is a circular 
spike at the vortex core; therefore, this contour plot is not shown. As the vortex interacts with the shock 
at T - 8, vorticity appears downstream. At T = 50, the vorticity becomes clearer. The vorticity patterns, 
like the convective density disturbances, look like the spokes of a wheel radiating from the vortex core 
and terminating at the shock. Looking closely at a region in the immediate vicinity of the vortex core at 
T = 50, it is clear that the lines of constant vorticity are oblong in shape (fig. 50). 

Velocity. Figures 51 and 52 show the changes in the axial and radial components of velocity. Except 
for small changes in amplitude, the velocity field shows no significant variation as a result of the inter- 
action on this scale. There are, however, features similar to those observed in the density perturbations 
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Figure 47. Mach-Zehnder interferogram of sound wave and contact surfaces generated shock and vortex interaction. (From 
ref. 46.) 



Figure 48. Contours of density perturbation downstream of shock at T - 50. 


when the contour range is decreased. Figure 53 shows the change in the axial component of velocity on 
a smaller contour range. The vortex, acoustic wave, and entropy disturbances are all readily distin- 
guished on this scale. 

Entropy . Figure 54 shows the change in entropy, which is defined as bs = s Q (p/p y )~ s 0 . At the initial 
time, there is no fluctuation in entropy, since the initial vortex is isentropic. Therefore, the contour plot 
at T = 0 is not shown. However, as the vortex core interacts with the shock wave, entropy appears 
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Figure 49. Contours of vorticity perturbation downstream of shock at T = 8 and T = 50. 
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Figure 50. Contours of vorticity perturbation in region immediately surrounding vortex filament at T= 50. 


downstream. Note that the entropy perturbations correspond to the convective fluctuations in density 
and vorticity as shown in figures 46 and 49. 

Flow features. Figures 40 through 54 illustrate flow features that result from the interaction of a 
ring vortex and a shock wave. After observing the perturbations in the flow quantities, it is clear that an 
acoustic wave and contact surfaces result from this interaction. The acoustic wave is isentropic and 
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Figure 51. Contours of axial velocity perturbation of shock at T = 0, T= 8, and T = 50. 

travels at the sound speed relative to the moving fluid. The contact surfaces do not support a pressure 
disturbance, but do support differences in flow velocity, density, entropy, and vorticity. The distur- 
bances are the result of the interaction. From a mathematical perspective, the shock wave is the feature 
that allows for the coupling of the linear modes (acoustic, entropy, and vorticity). Crocco’s theorem pro- 
vides the relationship between the flow vorticity and entropy. From this theorem it can be shown that 
vorticity can be generated in the presence of a curved shock. Because the evidence points to the impor- 
tance of the shock wave in the generation of the acoustic, vortical, and entropy disturbances, a closer 
observation of the shock dynamics will be made in the next section. 
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Shock dynamics. Figure 55 shows a space and time diagram of the shock displacement. For times 
r=0 through T = 6 the shock position remains nearly fixed as the vortex convects supersonically 
towards it. Once the vortex core begins to strongly interact with the shock at T = 6, the shock motion 
becomes significant. Figure 55 illustrates the nature of the shock displacement as a function ot time. A 
small disturbance on the shock is initiated at T ~ 6. The disturbance proceeds to travel both away from 
and towards the axis of symmetry. The speed of these waves is approximately the same in both direc- 
tions, but the disturbance reaches a local maximum away from the axis of symmetry and continues to 
grow towards the axis of symmetry. 
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Figure 53. Contours of axial velocity perturbation at T = 50 with range of contour levels reduced to show velocity features. 
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Figure 54. Contours of entropy perturbation at T = 8 and T = 50. 


Figure 56 clarifies the details of the shock displacement at various times. This figure shows the 
axial displacement of the shock as a function of radial position along the shock for T= 1, T = 6, T= 10, 
T = 20, T = 30, T = 40, and T= 50. At T= 10, the shock wave has maximum displacement magnitudes 
of 0.1 and 0.14 in the upstream and downstream directions, respectively. At this early time, the dis- 
placement of the shock is limited to a small region around r = 125, which is the radial position of the 
vortex filament. As time progresses, the maximum shock displacement occurs farther away from the 
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Figure 55. Shock displacement as a function of space and time. 


Time 

1 

6 



Figure 56. Shock displacement for T = 1,7 = 6, and T = 10 through T = 50. Positive displacement refers to downstream shock 
displacement. 


interaction point. The maximum upstream displacement occurs at a time of T ~ 10. However, the dis- 
placement in the downstream direction continues to increase, even after 50 periods. If the motion ot the 
maximum shock is considered to be displacements along the shock wave, it is apparent that these distur- 
bances travel away from the vortex and shock interaction point at a particular velocity. By careful exam- 
ination of the results presented in figure 55, this velocity is found to be » 0.54 core radii per period 
(which corresponds to a normalized velocity of * 0.68 core radii per period). This velocity is well below 
both the upstream and downstream sound speeds (1.0 and 1.15, respectively). 

The shock wave displacement continues to increase in the direction towards the center of the ring, 
while it maximizes and then decreases in the direction away from the axis of symmetry of the ring. 

Figures 55 and 56 show the asymmetry in the shock position relative to an axis passing through the 
vortex filament position at r - 125. The portion of the shock closest to the axis ot symmetry continues to 
move downstream over the time of the calculation. The maximum shock displacement is greater down- 
stream (Ax s = 0. 16) than upstream (Ax s = -0. 1). 
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Figure 57. Shock displacement and density and pressure perturbations for T = 50. Density and pressure perturbations were 
obtained slightly downstream of the shock (jc = 7.02). 



Figure 57 shows the displacement of the shock and the perturbation in density and pressure immedi- 
ately downstream of the shock at jc = 7.02 at T = 50. This figure shows that there is a correlation 
between the disturbances downstream of the shock and the shock displacement. The shock displacement 
leads the disturbances in pressure and density. The maximum and minimum of the shock displacement 
correspond to the large jumps in density at r - 100 and r = 150. The changes in slope of the 8x s and r 
curve that are located at r = 70, r = 1 20, and r = 1 30 also correspond to significant high frequency 
changes in density. This type of fluctuation in pressure occurs only at r= 70. Thus, by comparing this 
figure with the contour plots of figures 40 and 46, the spoke-like disturbances in density and vorticity 
can be related directly to the shock motion and curvature. 

Frequency analysis . To obtain information about the frequency content of the solution, the fast 
Fourier transform (FFT) is used to determine the time dependent information obtained from the numer- 
ical calculation to the frequency domain: 

P(t o) = f~ bp(t)e~ iat dt (105) 

• —oo 

where 8 p is the perturbation pressure defined as the difference between the pressure and the mean pres- 
sure (8 p = p - /?q), t is time, and co is the cyclic frequency. 

This integral is written in discrete form 


N-i 

P « At ^ p(nAt)e /amA/ (106) 

n= 0 

where At is the discrete time interval of the discrete pressure time history p(nAt ), and N is the number of 
points in the time history. Equation (106) is computed using the FFT. To efficiently use the FFT, N must 
be a power of 2. 

In this numerical simulation, solutions were saved at every period (defined as T - rju ^ ). Thus, there 
are 50 discrete representations of pressure in the time history that are available for analysis at each point 
in the computational domain. To effectively use the FFT approach, the time history at each point at 
which the spectrum is to be computed is padded with zeros so that the number of discrete values in the 
time history is a power of two. Padding the time history with zeros does not affect the highest frequency 
that can be resolved in a discrete approximation to the Fourier transform (the Nyquist frequency), since 
the Nyquist frequency co c is a function only of the temporal increment in the time history (O c = n/At. 
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However, padding the time history with zeros does increase the frequency resolution of the spectral 
estimate. 


The increase in frequency resolution provided by zero padding is readily illustrated by an example 
(ref. 34). Suppose the original time history T has b elements: T(nAt), n = 0,1,2,...,/? - 1. The original 
time history is then padded with additional b elements, which are zero: T(nAt) = 0 ,n = b,b+ 1,...,2 b - 1 . 
The discrete Fourier transform T r is then 
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Note that = nk/bAt , k = 0,1,2,...,/?, and the frequency resolution is Ato = lUbAt. Similarly, for the orig- 
inal time history without zero padding, <£> k = InkJbAt, k = 0,1,2,..., b/2, and the frequency resolution is 
Ato = 2n/bAt. Thus, the zero padding has the effect of increasing the frequency resolution by a factor of 
two. 


Although it would be sufficient to pad the 50-element time history with 14 zeros to obtain a 64 ele- 
ment time history, each time series was padded with 206 zeros for a total of 256 points in the time his- 
tory. This number of zeros provided better frequency resolution that was necessary for localization of 
the acoustic energy in the spectrum. The price paid for enhanced frequency resolution is a loss of accu- 
racy. (See ref. 34 for details.) 

Typical sound pressure levels arc presented in figures 58 and 59. The frequencies are normalized by 
the upstream flow velocity and core radius to obtain a Strouhal number. In these figures, the effect of 
the mean flow has been removed from the spectra by subtracting the time average of the pressure from 
each point in the time history. This subtraction removes the energy from the zero frequency bin. The 
sound pressure level is presented in decibels (dB). The decibel level is obtained by computing 


P( to) = 20 log 


- Pref - 


( 108 ) 


where p K f = 20 [iPa is the conventional acoustic reference pressure. 



Figure 58. The SPL at 6.081 1 core radii from source and 45° from horizontal as a function of dimensionless frequency. 
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Figure 59. The SPL at 1 1 .74 core radii from source and 45° from horizontal as a function of dimensionless frequency. 



Figure 60. The SPL at 23. 1 core radii from source and 45° from horizontal as a function of dimensionless frequency. 

Figures 58 through 60 show the sound pressure level at locations along a 45° line from the sound 
source point (x = 7, y = 125) at 6.08, 1 1.74, and 23.1 core radii away from the source point. The highest 
energy levels are in the low frequency range (Strouhal numbers less than 0.1). There are additional 
peaks at Strouhal numbers of = 0.2, 0.4, and 0.6. The energy in the frequency band of about 0.22 
becomes relatively more important than the low frequency energy as distance from the source point is 
increased. By observation of the time history from which the spectrum is computed, the period of the 
acoustic wave is determined to be approximately 21 periods. This period corresponds to a Strouhal 
number of to r c /u x = 2nr c /Tu x = 2711/21(1.5) = 0.2. Thus, the energy in the frequency range of 0.22 is 
associated with the acoustic energy of the signal. The large peaks near the Nyquist frequency (Strouhal 
number = Jt) are probably due to aliasing. Hardin has shown that for spectral estimates near the Nyquist 
frequency co c contributions from -(D c can appear, even for a sufficiently sampled time history (ref. 57). 
Figure 61 shows the sound pressure level of the two frequencies 0.025 and 0.2 as a function of distance 
from the sound source point. The data are taken along a 45° line from the horizontal that passes through 
the vortex core at r = 125. The data show a rapid decay in the energy associated with the low frequency 
energy. The energy decays at about 6 dB per doubling of distance away from the source. This implies a 
decay rate of 1/r 2 , which is consistent with the expected decay rate for a vortical pressure field. 

The energy at the Strouhal number of 0.22 decays at a lower rate of = 3 dB per doubling of distance, 
which implies a sound pressure level decay of 1/r. Because the sound pressure level is a function of the 
square of the pressure, this result is consistent with the pressure decay rate of l/r 1/2 that was discussed 
in the section entitled “Pressure” and implies cylindrical spreading of the acoustic energy. 
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Figure 61. The SPL at 45° from point where a horizontal line passing through vortex filament (r = 125) passes through undis- 
turbed shock. 


Sound intensity level . The sound intensity is a vector quantity defined by 

I(X) = f f w dt (109) 

1 Jo 

where T is the period of a periodic disturbance or a sufficiently long time for nonperiodic bounded sig- 
nals and W is the instantaneous acoustic energy flux vector 

W = (5p + p 0 U 0 -5u)f8u + ^uJ 
v p 0 J 

where U is velocity and 5 p, 5 w, and 5p are disturbances in pressure, velocity, and density, which are 
defined 


bp = p-p 0 

5u = u - u 0 
Sp = P-P 0 

and the subscript 0 represents the mean flow state. 

Note that in a quiescent medium (where U 0 = 0), three of the terms in the instantaneous acoustic 
energy flux vector vanish and the equation for sound intensity reduces to 

I c (x) = ^ | bpbu dt (110) 

Because early acoustics theory was developed for sound disturbances in a quiescent state, this defi- 
nition will be called the “classical” definition of sound intensity and it will be labeled with a subscript c 
to distinguish it from the definition of intensity presented in equation (109). It will be shown later that it 
can be beneficial to consider the classical definition of sound intensity, even in the presence of mean 
flow. 

In the calculations performed here, the acoustic signal is transient, which provides some difficulty 
in the interpretation of sound intensity. Because this problem models the periodic convection of turbu- 
lent disturbances through a shock wave, the interpretation taken here is to use the length of time of the 
transient disturbance as the period of a signal. This interpretation assumes that sound waves generated 
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by periodically convecting disturbances do not overlap and the time between pressure disturbances is 
equal to the time of the transient. 

Once the sound intensity is computed and dimens ionalized, the sound intensity level is determined 
by the relation 


IL = 10 logy— 

'ref 

where / ref = 10 watts/m . To dimensionalize the calculations, standard sea-level conditions are used; 
the upstream stagnation sound speed is 340 m/s and upstream static pressure is 101 325 N/m 2 . 

It is important to state that the computation of the sound intensity vector will not guarantee that the 
result is all sound. It is difficult, if not impossible, to separate the acoustic fluctuations from hydrody- 
namic fluctuations in the near field of the sound source. The definition of the sound intensity vector can 
provide useful information regarding the strength and directivity of disturbances in the flow and will 
provide a useful description of the sound when the calculation is made in the far field. 

Typical intensity plots are presented in figure 62, which shows the axial and radial components of 
the intensity vector for an upstream Mach number of 1.5 and a vortex strength of 0.75. The total time 
history in the calculation of the intensity components is T = 50. 

These plots show two distinct regions of high intensity level. One region is along the path that the 
vortex travels after intersecting with the shock wave. This region is most clear on the plot of the axial 
component of intensity level. This region is not sound intensity, since most of the disturbances related to 
the vortex convect. The high intensity level shown in this region demonstrates the high correlation 
between the pressure and velocity disturbances in the vortex. 

Some acoustical energy may be generated as the vortex changes shape after its interaction with the 
shock, but small scale vortical motions were not observed to be a significant sound source in these cal- 
culations. This insignificance is most clearly seen in figure 43, which shows a carpet plot of pressure 
fluctuations downstream of the shock. The most significant structure is the ring, which has already been 
identified with the interaction of vortex core with the shock wave. Much less significant pressure distur- 
bances are seen in the center of this acoustic wave, but there are no waves visibly emanating from the 
vortex. 

The other region of high intensity is visible along the shock wave. The intensity plots show that 
sound is directed primarily in directions closely aligned with the shock wave. The intensity originates at 
the point of interaction and travels both towards and away from the axis of symmetry. Both axial and 
radial components of intensity have high levels along the shock; the region is much narrower in the con- 
tour plot of the axial component. The high amplitude of these waves makes it difficult to classify them 
as acoustic in the classical sense, but because there is no convective velocity in the radial direction, 
these disturbances are clearly not convective. 

A plot of the axial component of sound intensity that is determined by the classical definition is 
shown in figure 63. (Because the mean flow is in the axial direction only, the radial component of clas- 
sical sound intensity is the same as that shown in fig. 62.) The region of high intensity in the vicinity of 
r = 125 is due to the passage of the vortex along this path. Figure 63 shows that the intensity has four 
lobes not related to the convection of the vortex. Two of these lobes are along the shock wave, similar to 
the results shown in figure 62. However, two additional lobes originate at the point of vortex filament 
and shock interaction and point at angles of approximately 50° and -55°. The lobe directed at —55° is 
narrower than the lobe directed at 50°. Note that these regions are not evident in the contours of the full 
description of the sound intensity (fig. 62) because the energy associated with the acoustics is over- 
whelmed by the contributions from the mean flow terms. The high intensity level along the lobes is 
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Figure 62. Sound intensity level. 


significant because it corresponds to the location ot the strongest region of the sound wave. (See 
fig. 42.) This correspondence indicates that it can be beneficial to consider the classical definition of 
sound intensity, even in the presence of a mean flow. It is legitimate to apply equation (1 10) as long as 
it is not used to draw conclusions about acoustic energy conservation. 

Effect of Mach number on directivity. To determine the effect of Mach number on the directivity of 
the pressure disturbances, a series of computations was performed where all flow parameters were held 
fixed except the upstream Mach number. The preshock Mach numbers studied were M - 1.1, 1.3, 1.4, 
1.5, and 1.7. These Mach numbers were chosen because they are within the range of practical interest. 
The sound intensity level was computed over this Mach number range. The directivity angles as a func- 
tion of shock strength are plotted in figure 64. The figure shows the angles of the four primary beaming 
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Figure 63. Sound intensity level with primary directivity of sound wave along shock wave and downstream at angles of ±50°. 
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Figure 64. Directivity angles of primary intensity lobes as a function of upstream Mach number. 


directions. The change in lobe direction is most significant for the lobes closely aligned along the shock. 
As the Mach number increases, these lobe angles approach ±90°. 


Effect of flow Mach number on sound pressure level Because of the complexities in defining and 
computing the sound intensity level for a transient signal in the presence of a mean flow and because the 
human ear responds to pressure fluctuations, a study of the effect of flow Mach number on the sound 
pressure level is presented here. In this study, sound pressure level is computed for a point in the acous- 
tic far field (defined to be one wavelength away from the source) at a distance of 23 core radii from the 
interaction point at an angle of 0 = 50°. In these calculations, the sound pressure level is computed by 


SPL = 20 


( 


log 


^rms 

\Pre( j 


(HI) 
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where p ref is standard reference pressure and 



where T is the period of the acoustic signal. Because this calculation is performed in the far field, the 
acoustic wave is readily distinguished from hydrodynamic disturbances. The result of this study is pre- 
sented in figure 65, which shows sound pressure level as a function of P = JM 2 - 1 , for Mach num- 
bers in the range 1.1 < M < 1.5. Experimental measurements obtained by Seiner and Norum for shock 
noise of an underexpanded supersonic jet, which was measured at a distance of 12 ft from the jet center- 
line (145.6 jet radii from the source), but corrected to a distance of 0.188 jet radii for comparison with 
numerical computation are also presented in figure 65 (ref. 58). In addition, the trend that the sound 
pressure level of shock noise in imperfectly expanded supersonic jets is proportional to p 4 over a large 
range of flow Mach numbers, which was observed by Harper-Bourne and Fisher is also included in the 
figure (ref. 63). The slope of the sound pressure level and p 4 curve very nearly matches the slope of the 
experimental results for shock noise in supersonic jets. This similarity is significant because it shows 
that this simple model for shock noise generation can predict the effect of Mach number on shock noise 
that is observed in experiment. 

Although this model correctly reproduces the trend in sound intensity level as a function of shock 
strength, it does not reproduce the actual sound amplitude even when differences in the distance 
between the sound source and measurement location are accounted for. This lack of reproduction is not 
surprising since shock noise measured during an experiment is the result of many interactions of turbu- 
lent structures of a variety of sizes and strengths with a sequence of shock waves of decaying strength 
and at varying angles to the flow. 



Figure 65. The SPL as a function of upstream Mach number. 


Strong Interaction 

The results presented in this section are for the interaction of a strong vortex with a shock wave. For 
all the results presented in this section, the upstream Mach number is 1 .5 and the strength of the vortex 
is T = 5.5. Thus, this vortex is 7 1/3 times stronger than the vortex studied in the previous sections. 
When the vortex strength increases, the pressure disturbances generated downstream of the shock 
increase in magnitude, which shows contour plots of the pressure disturbances at T = 50 (fig. 66). The 
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Figure 66. Contours of pressure perturbation downstream of shock at T = 50 and V = 5.5. 

actual range of the pressure disturbances are on the order of 1, but for the plot, the range of contours was 
chosen to be from -0.01 to 0.007 so that direct comparison with figure 40 can be made. In both figures, 
a sound wave and a cusp-shaped pressure disturbance are apparent on the side of the vortex filament 
closest to the axis of symmetry and downstream of the shock. However, the sound wave is much stron- 
ger in the strong vortex and shock interaction case, and on this plot contour scale, additional pressure 
structures are visible downstream and below the vortex filament (fig. 66). 

Figure 67 shows pressure perturbation along radii extending from the vortex center at 10° incre- 
ments. This figure shows that the peak-to-peak pressure amplitude is a maximum at 55° and -55°, as in 
the case for V - 0.75. From this figure, the peak- to- valley measure of the wavelength is found to be 
= 2.5. In comparison with figure 42, the disturbances tor the strong vortex interaction case are almost an 
order of magnitude larger. For example, the amplitude of the peak disturbance at 50° is 10.2 times larger 
for the strong interaction case (T = 5.5) than the weak interaction case (T = 0.75). Thus, the acoustic 
pressure p scales as the vortex strength T. 

Note that the pressure disturbances immediately downstream of the shock form very steep gradi- 
ents. These gradients may be considered to be shock waves and have been referred to as “reflected 
shocks” (ref. 10). Figure 68 shows the distribution of perturbation pressure as a function of distance 
from the initial shock position at locations above and below the vortex filament. There are significant 
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Figure 67. Pressure perturbations along radii extending from vortex core with f = 5.5 at a = 30. 


74 




Figure 68. Pressure perturbations as functions of axial position at T- 50. 


jumps in pressure downstream of the original shock wave, especially below the vortex filament. The 
formation of steep gradients downstream requires that the algorithm used in the calculation be robust. 

Figure 69 shows contours of the density perturbations downstream of the shock. The same features 
present in the T = 0.75 case are present here, but the strength of these features is now enhanced. 

Typical Interaction of Clockwise Vortex 

In this section, results are presented for the interaction of a clockwise rotating vortex with a shock. 
This case is not representative of the physics of the interaction of vortices in the jet shear layer with 
shock waves, but is included for completeness. This interaction could model wake flow and shock inter- 
action, such as the wake of a helicopter blade interacting with a shock on the subsequent blade. 

Figure 70 shows contours of pressure for the vortex rotating clockwise with a shock wave. The 
strength of the vortex is V = -0.75 and the Mach number upstream of the shock is M = 1 .5. 

Note that this case is directly analogous to the case presented in the section entitled “Strong Interac- 
tion” except for the sign of the vortex circulation. Comparison of figure 70 with figure 40 shows that the 
contours of pressure perturbation are quite similar except that the regions of compression in figure 70 
are regions of rarefaction in figure 40. The difference in the rotation sense of the vortex results in a dif- 
ferent response of the shock wave, which in turn results in a difference in the sign of the pressure distur- 
bance downstream. As shown in figure 55, a CCW rotating vortex causes the shock to move upstream 
in the region above the votex filament and to move downstream in the region below the vortex filament. 



Figure 69. Contours of density perturbation downstream of shock at 7= 50 and f = 5.5. 
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Figure 70. Contours of pressure perturbation downstream of shock at T - 50. 


For the clockwise rotation, the shock wave moves downstream in the region above the vortex filament 
and upstream in the region below the vortex filament, as shown in figure 71. Figure 72 clearly shows 
that the shock displacement is significantly greater upstream than downstream. 

The downstream pressure disturbances resulting from the CW and CCW vortices are not perfect 
images of each other about the vortex filament. The cusp-shaped structure described in the section 
entitled “Strong Interaction” connects the acoustic disturbance to a position on the shock near the axis 
of symmetry in both cases. 

There is also asymmetry of the acoustic disturbances. Figure 73 shows the pressure perturbations 
along radii extending from the vortex core at ±40°, ±50°, and ±60°. For ease of comparison, the results 
have been plotted so that the pressure disturbance associated with the positive angle of the CW vortex is 
compared with the negative angle of the CCW rotating vortex. The results show good agreement 
between the pressure peak at the positive angle for the CW vortex and the negative angle for the CCW 
vortex. However, there is a significant difference in the maximum pressure amplitude of the distur- 
bances that corresponds to the CW vortex at a negative angle and the CCW vortex at the positive angle. 
For these cases, the CW vortex generates a larger pressure disturbance. The largest difference is 
156 percent for the sound disturbance traveling at ±60° from the horizontal. 
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Figure 71. Shock displacement as a function of space and time. 
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Figure 72. Shock displacement as a function of radial distance for 7 = 1,7=6, and 7=10 through 7 = 50. 
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Figure 73. Pressure perturbations along radii at ±40°, ±50°, and ±60° extending from vortex core at jc = 30, y - 125. 




Effect of Vortex Core Size 

To show the effect of vortex core size relative to the vortex ring size, a study is included in which 
the ratio of the core radius to ring radius is 1/250. Figure 74 shows the contours of change in pressure at 
T = 50 for this case. The range of pressure levels in the contour plot is kept the same as in figure 40 for 
direct comparison. The strength and wavelength of the acoustic waves are essentially identical to those 
in the case where the ratio of core radius to ring radius is 1/125. The primary difference is that the cusp- 
shaped structure so apparent in the core radius to ring radius shown in figure 40 is no longer visible on 
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Figure 74. Contours of pressure perturbation for r=0.75, upstream M = 1.5, and ratio of vortex core radius to ring radius is 

1/250. 


the contour plot. This lack of a cusp-shaped structure indicates that the presence of this wave may be 
related to axisymmetry. As the vortex ring becomes larger relative to the core size, the resultant wave 
structures more closely resemble those generated in a two-dimensional planar interaction, which does 
not show the cusp-shaped structure. 

Summary of Chapter 4 

In this chapter, the interaction of a ring vortex with a shock wave is presented as a simple model for 
a mechanism responsible for generating shock noise in imperfectly expanded supersonic jets. This 
model was inspired by the works of Ribner (refs. 49 and 50) and Moore (ref. 38) who considered two- 
dimensional interactions of disturbances with shock waves in the context of linear theory. The model 
reported herein differs from these pioneering studies because the interaction considered more closely 
models the interaction of axisymmetric turbulent structures within axisymmetric jets. In addition, the 
present study provides a tool tor solving the nonlinear equations that govern the fluid dynamics of this 
very complex interaction. 

An effort was made in the design of these calculations to model the physical parameters of an 
imperfectly expanded jet. The size of the vortex core relative to the vortex ring radius was small because 
of numerical considerations, so flow parameters were chosen to model the interaction of a turbulent dis- 
turbance with a shock wave in the first shock cell. The magnitude of the vortex strength and the sense of 
the vortex rotation were chosen to be consistent with those observed in experiments. The range of Mach 
numbers studied was within the range for practical nozzles. 

Observation leads to the conclusion that sound, entropy, and vorticity are generated downstream of 
the shock wave during the interaction process. The sound wave is isentropic and propagates down- 
stream at the sum of the convection and sound speeds. Contact surfaces are formed during the interac- 
tion process. These disturbances do not support a difference in pressure, but do support differences in 
density, entropy, and vorticity. Both sound waves and contact surfaces have been observed in experi- 
mental studies of shock and vortex interaction. It is believed that these are the first calculations per- 
formed with enough resolution to show the presence of the contact surfaces downstream of the shock 
for these flows. Contact surfaces have been observed in experiments. These experiments provide evi- 
dence for the physical nature of these disturbances, but as discussed in chapter 3, numerical error asso- 
ciated with moving shock waves may readily manifest itself in entropy, so further analysis is necessary 
to validate the strength of the computed contact disturbances. 

The structure of the contact surfaces is related to the shock dynamics. The contact disturbances are 
observed to be generated by a wave that is initiated during the interaction and travels along the shock 
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both towards and away from the center of the vortex ring. The motion of the shock corresponding to the 
wave traveling along its length generates entropy and vorticity disturbances. Pressure perturbations are 
also observed to originate at the locations of large shock motion. 

Analysis of the data provides further insight into the physics of the interaction. A frequency analysis 
suggests that the sound wave generated in these computations decays as a cylindrical wave. This cylin- 
drical decay rate is due to the presence of the shock wave, which acts as a barrier to sound spreading 
upstream of the shock. An analysis of sound intensity level provides insights into the directive nature of 
the sound wave. High intensity levels are observed close to the shock wave and in downstream direc- 
tions at angles of approximately 50° and -55° from a horizontal axis through the vortex filament for the 
case where the upstream Mach number is 1.5. 

Additional studies of the interaction of a strong vortex with a shock wave validated the robustness 
of the numerical code used in this work. The interaction produced extremely high gradients in the flow 
downstream of the shock and the computation maintained stability. 

The sound generated by clockwise and counterclockwise rotating vortices was compared. It was 
found that the structure of the alternating compression and rarefaction zones along the sound wave 
changed sign when the sense of the vortex rotation changed. (Regions of compression resulting from the 
clockwise vortex became regions of rarefaction for the counterclockwise vortex and vice versa.) It is 
also observed that the clockwise rotating vortex generates disturbances that can be significantly larger 
in amplitude. The peak pressure perturbation resulting from the interaction of the clockwise rotating 
vortex and the shock produces pressure levels up to 156 percent higher than the counterclockwise rotat- 
ing vortex. 

A study of the effect of the ratio of core size to ring size was also performed. This study shows that 
the effects of axisymmetry are reduced when the core radius size is decreased relative to ring radius. 

A study of the effect of flow Mach number on the sound pressure level found that the sound pres- 
sure level increased with Mach number. The rate of this increase corresponds closely to that observed in 
experimental data of shock noise of supersonic jets. This result implies that the interaction of a vortex 
ring with a shock wave is a dominant physical process in shock noise generation in supersonic jets. 

Conclusions 

This research pioneered the application of a direct computational approach to the study of shock 
noise mechanisms. Direct simulation of sound generation in shocked flows is challenging because of the 
disparity in amplitudes between the acoustic waves and the shock waves. These challenges were met by 
the implementation of an essentially nonoscillatory (ENO) scheme that uses adaptive stenciling to 
maintain high-order accuracy in smooth regions of the flow. This maintenance of high-order accuracy 
minimized numerical dissipation of the acoustic waves while maintaining sufficient numerical dissipa- 
tion at the shock for stability. A study of the economics of high-order schemes shows that the added cost 
of higher order algorithms should be balanced with the level of accuracy required. An analysis per- 
formed for sound in a converging nozzle showed that for a numerical error on the order of 1 (T 6 , the 
third order accurate ENO scheme is the most cost-effective. Therefore, a third order ENO algorithm was 
used for most of the work presented herein. A study of the numerical error in the computation of slowly 
moving shock waves showed that spurious numerical waves were produced downstream of the shock. 
The numerical error manifested itself primarily in entropy and is a function of the algorithm used in the 
computation and the shock speed relative to the grid. As the shock speed relative to the grid increases, 
the entropy error decreases. 

This report presents and describes the modeling of sound generating mechanisms in a supersonic 
jet. Experimental evidence is presented to illustrate that shock noise contributes significantly to the 
sound field of a supersonic jet. Two mechanisms of sound generation by shocked flows were 
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investigated: shock motion and shock deformation. These mechanisms were modeled by the interaction 
of sound disturbances with shock waves and the interaction of vortical structures with shock waves. 
These models permit the consideration of shock oscillation and shock deformation in the development 
of sound. 

Analysis of shock motion with Lighthill’s equation showed that monopole, dipole, and quadrupole 
terms all have potential to contribute to the far field sound. At low supersonic Mach numbers, the mono- 
pole term dominates, followed by the dipole. The dipole term is highly directional. 

Shock motion was modeled numerically by the interaction of a sound wave with a shock wave. Dur- 
ing the interaction, the shock wave begins to move and the sound pressure is amplified as the sound 
wave passes through the shock wave. Computations of sound waves interacting with shocks in a 
converging-diverging nozzle were performed. The results show that the amplitude of the transmitted 
pressure perturbation is greater than the incident pressure perturbation for all Mach numbers greater 
than one. The numerical approach was validated by comparison of the computed ratio of transmitted 
with incident sound pressures with linear theory. The comparison is good over the range of shock Mach 
numbers studied (1.3 < M < 2.6). An energy analysis was performed to determine whether acoustic 
energy was generated in the sound and shock interaction. The analysis is based on an exact representa- 
tion of the transport of energy in an arbitrary flow field and shows that acoustic energy is generated dur- 
ing the sound and shock interaction. The source term of this energy was shown to be confined to a 
region along the shock wave in space and time, and is a function of the changes in entropy, momentum, 
and temperature in the mean flow state. 

Shock deformation was investigated by the simulation of a ring vortex interacting with a shock 
wave. This investigation has practical significance because it models the passage of a turbulent structure 
through a shock wave. Observations of the evolution of perturbations in pressure, density, entropy, vor- 
ticity, and velocity downstream of the shock lead to the conclusion that acoustic waves and contact 
surfaces are generated by the interaction of shock and vortex. That these two fluid structures were gen- 
erated by the interaction was validated by experimental evidence. The structure of the contact surfaces 
was related to the shock dynamics. The contact surfaces were observed to be generated by disturbances 
that were initiated during the interaction and traveled along the shock both towards and away from the 
center of the vortex ring. The motion and deformation of the shock generated entropy and vorticity, 
respectively. 

Analysis of the numerical results demonstrated that the sound wave that results from the interaction 
of a vortex ring with a shock wave spreads cylindrically. This cylindrical spreading is due primarily to 
the presence of the shock wave, which acted as a barrier to sound traveling upstream. Analysis of the 
sound intensity level over the region of the computation provided insight into the directivity of the 
sound. High intensity levels were seen along the shock wave at angles of approximately 50° and -55° 
from a horizontal axis through the vortex filament when the upstream Mach number was 1 .5. The peak 
sound amplitude that was generated by a clockwise rotating vortex was found to be as much as 
156 percent higher than sound generated by the interaction of a counterclockwise rotating vortex and 
shock wave. A significant result of this work is that the sound pressure level was shown to increase with 
shock strength. The relationship between the sound pressure and shock strength, defined by the parame- 
ter J3 = 1 , is shown to be approximately: sound pressure «: p 4 . This relationship is consistent 

with experimental observations of shock noise in supersonic jets and implies that the interaction of a 
vortex ring with a shock wave is a dominant physical process in the physics of shock noise generation. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
November 8, 1996 


80 



Appendix A 

Derivation of Unsteady Shock Jump Relations 

In this appendix, the shock jump relations for a moving shock are derived by using generalized 
functions. The results obtained here for the continuity and momentum equations are also in reference 9 
without a derivation of the unsteady shock jump relation for the energy. 

Generalized functions will not be defined in a rigorous mathematical context here, but will be pre- 
sented to show the connection between generalized and ordinary functions for clarity. (See ref. 9 and its 
references for details.) 

Conventionally, a function is defined as a table of ordered pairs (x, f(x)), where for each x, f(x) is 
unique. This table may have an infinite number of ordered pairs. In an analogous fashion, in generalized 
function theory, the function /(x) is defined by its action on a given space of ordinary functions called 
test function space: 


FI*] = /(x) <t>(x) dx (Al) 

J — oo 

where the function ty(x) is a test function that must satisfy certain properties given in reference 9. The 
mapping described by equation (Al) is functional. The function /is now identified by the new table 
F[<|>], <|> g test function space. 

It can be shown that there are an infinite number of functions <j> that satisfy the conditions prescribed 
on the space of test functions, so the table produced by equation (Al) has an uncountably infinite num- 
ber of elements. The space of test functions is so large that the functions on this space generated by 
equation (Al) contain not only ordinary functions, but also generalized functions. Thus, ordinary func- 
tions are a subset of the generalized functions. It can be shown from classical Lebesque integration the- 
ory that the Dirac delta function cannot be an ordinary function (ref. 9). However, functions such as the 
Dirac function are included in the definition of generalized functions. Now that the concept of function 
space has been extended to include functions such as the Dirac function, the extension of the definition 
of derivative is presented. 

Let f(x<,y,z) be a piecewise smooth function with one surface of discontinuity. Denote this surface by 
8 = 0. Such a surface is illustrated in figure Al. At this surface of discontinuity, there is a jump in the 
value of the function denoted by 


4/ = /(g = 0V/(g=0 ) 


(A2) 



Figure Al. Schematic of a discontinuous surface. 
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Note that g = 0 + is on the side of the surface into which Vg points and that, following Farassat, S7g = n 
(ref. 9). 

The generalized divergence of a vector /V ■ f is related to the ordinary divergence V f and the 
jump across the surface normal Af by 

V f = V • f + Vg • Af5(g) 

= V • f + n • Af8(g) (A3) 

The shock jump relations for a flow in which a discontinuity in the flow is moving may be derived 
by using equation (A3). The differential forms of the laws governing the conservation of mass, momen- 
tum, and energy are valid even when discontinuities exist within the region, as long as the derivatives 
are interpreted as generalized derivatives. The application of generalized derivatives to the conservation 
laws is most readily performed when the equations are in divergence form. Begin by considering the 
conservation of mass 


^ + V-(pu) = 0 (A4) 

Applying the definition of the generalized derivative equation (A3) 

3* + A(p >§7 + V ' (pu ^ + A(pu ^ ‘ n5(g ) = 0 (A5) 

The sum of the first and third terms defines the ordinary continuity equation, thus this sum is zero. 
The term dg/dt is equal to the negative of the velocity of the surface — v n . 

Thus, 


A[p(w n -v n )] = 0 (A6) 

where u n = u n is the component of velocity normal to the surface g and v n is the velocity of the sur- 
face. If the states upstream and downstream of a shock are denoted by the subscripts 1 and 2, respec- 
tively, and the shock is not moving, note that equation (A6) reduces to the Rankine-Hugoniot relation 
for steady flows in one dimension: 


p 2 u 2 “ P i w i 

Consider now the conservation of momentum: 


d(P»,) jkp «f«/) dp _ 

dt + dxj dx ; 


Applying the definition of generalized derivatives 


3(Pm) d(p «,«■) 

— + I— £- H — + 

dt dxj dx t 




8(f) = 0 


(A7) 


(A8) 


But, since dg/dt = —v n and dg/dxj = tij, and the sum of the first three terms is equivalent to the 
momentum equation for continuous flows, equation (A8) reduces to 


Mpu i (u n -v n ) + pn t ] = 0 


(A9) 
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Again, note that when the shock is not moving, this equation reduces to the Rankine-Hugoniot rela- 
tion for steady, one-dimensional flow: 


2 2 
p 2 « 2 + />2 = Pi W 1 + P\ 


Now consider the energy equation 


de d(e + p)u i 

di + a*,. =0 (A10) 

where e is the total energy defined by e = ph + 1/2 p u 2 -p = pH-p. The total specific enthalpy is rep- 
resented by the symbol H. Applying the definition of the generalized derivative 


| + ^ + { 4w | +4[(e + pK ,i} 5(s) = 0 


Making simplifications similar to those made for the continuity and momentum equations 

-A(e)v n +A[(e + p)u n ] = 0 

-A(p/i + l -pu 2 - p + A ^p/i + X -pu 2 = 0 


A li ph + \ P “ 2 ) + PV n = 0 


For steady one-dimensional flow, equation (A 14) reduces to the Rankine-Hugoniot relation 


, 1 2 , 1 2 
h 2 + -p 2 u 2 = h x + ^9\ u \ 
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Appendix B 

Equations for the Velocity and Pressure of a Ring Vortex 

For the purpose of prescribing the initial condition for the numerical calculation, the vortex ring is 
assumed to be incompressible. The vortex moves relative to a fixed coordinate system at a velocity 
equal to the sum of the mean flow velocity U and the vortex translational velocity V. A cross section of 
the vortex ring is illustrated in figure Bl. The equations for velocity and pressure in the fluid as a result 
of the presence of the vortex are derived for the purpose of prescribing an initial condition for a numer- 
ical calculation. 


Velocity Outside the Core 

Lamb provides the expression for the stream function \\f of a vortex ring (ref. 59) 


V = - 



(Bl) 


r 



Figure B 1. Ring vortex moving at a velocity U+V with respect to a Fixed coordinate system (jc, r). 
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where T is the circulation, r is the distance from the axis of symmetry, and K and E are the complete 
elliptic integrals of the first and second kind 


where m is an index and 


2 2 
r\ \ 


K(m) - P (1 -m sin 2 9) 

J 0 


dQ 


E(m) = p ( 1 - m sin 2 9) 

J 0 


d 9 


k 2 = 


4rr„ 


[(*-* 0 ) 2+ ( r+r on 


(B2) 


(B3) 


(B4) 


where and xq are the radial and axial positions of the vortex filament, respectively. 

Given the stream function, axial and radial velocity components u and v, respectively, can be found: 


u - 


1^ 

r dr 


r 

2nr 2 


K(k) + 


(*-*o) ~ r 


E(k) 


(B5) 


ld\j/ 
r dx 


r 

27ir 2 rr 2 


K(k)- 


’ 2 , .2 2 ' 

r 0 + (x-x 0 ) +r 


E(k) 


where 


(B6) 


r \ = V( r - r o ) 2 (^ — ^o ) 2 

I 2 2 

r 2 = V^ r + r o) 

In the implementation of these equations for vortex velocity, the polynomial approximations of K(k) 
and E(k) y which have an error bounded by 2 x 1CT 8 , are used (ref. 64). 


Velocity Inside the Core 

To avoid mathematical singularity on the vortex filament and to better model the physics of a real, 
viscous vortex, the velocity distribution inside the vortex core is assumed to have a linear distribution of 
tangential velocity: 


u = 



(B7) 
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(B8) 


where v e is the tangential velocity at the core radius r c and R and 0 are the polar coordinates centered on 
the vortex filament. The tangential velocity v 0 is known by substituting x = x 0 + r c where r = r 0 into 
equation (B6). Note that this description of the core assumes that the core is circular. This equation is an 
approximation since the core is elliptical for finite values of r c lr$. 

Pressure Outside the Core 

The pressure field of this vortex is determined by momentum conservation. Consider the momen- 
tum equation for an irrotational flow: 





(B9) 


where p is density, p is pressure, q is the magnitude of the velocity vector, and u is the velocity vector 
[u,v] T . Now, the velocity component in the axial direction is u = U + u v (x - 2;(f),r) where U is the mean 
flow velocity, u v is the perturbation velocity induced by the vortex, the coordinate £, is the vortex fila- 
ment position, t is time, and r is radial distance. The radial component of velocity is v = v v (x-%(t),r), 
where v v is the velocity induced by the vortex. 

The flow is not steady in the fixed reference frame because the position of the vortex varies as a 
function of time h, = x + (U 4 - V)t, where V is the vortex translational velocity. Evaluating the derivative 
du/dt results in 


du 

dt 


du v (x-Z,(t),y) 

dt 

d[*-^(f)] dt 
du 

du v 

-97 (t/ + V) 


Similarly, 


di 




Note that for irrotational flow, dv/dx = du/dr. Therefore, dv v /dx = du v /dr so that 


(BIO) 


(B 11) 


du 

dt 


= -V uJU + V) 


(B12) 


If density is assumed to be constant, the momentum equation can be written 

2 


V [p + 2‘- (t/ + V)M ‘] = ° 


(B13) 
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where q 2 = u 2 + v 2 = ([/ + w r ) 2 + v 2 . Substituting for w and v and evaluating the integral of equation 
(B13) along a streamline that terminates in the mean flow state where the velocity induced by the vortex 
is zero, results in 


P~Po 


P 0 , 2 2, 

-yK + + Po Vw v 


where the subscript 0 denotes the mean flow state. 


(B14) 


Pressure Inside the Core 

Equation (B14) is valid for the flow outside the vortex core, where the flow is irrotational. To deter- 
mine the pressure field inside the vortex core where the flow is rotational, consider the radial momen- 
tum equation 


2 

pv fi 

dp — — - dr (B15) 

r 

Substituting for v$ and integrating from the edge of the vortex core to r, 

P(r ) = P(r c ) + ~\[r 2 - rfj (B16) 

r 2 

assuming that p = p 0 . The pressure p(r c ) that is determined from equation (B14) to ensure continuity of 
pressure across the vortex core interface is 

Pn 2 

p( r c ) = Po~y Vq (B17) 

The pressure profile of a counterclockwise vortex of strength T = 0.75 is shown in figure B2. Note 
the asymmetry in the pressure field relative to the vortex filament and the region of slightly positive 
pressure below the filament. 

When the vortex ring is the initial condition, the flow quantities are normalized with respect to the 
upstream static pressure and density and the vortex core radius. 

Remarks 

Note that although the pressure and the velocity are defined to be continuous at the rim of the vor- 
tex core, the derivatives of these quantities are not. Because the natural dissipation in the algorithm will 
automatically smooth these derivatives at the rim, discontinuities in derivatives at the rim are not a 
problem with this algorithm. Less robust schemes may have difficulty with this initial condition. 
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Appendix C 


Finite Wave Boundary Conditions 


General Description 

This report uses the inflow and outflow boundary conditions developed by Atkins and Casper 
(ref. 44). Some equations from reference 44 are examined in this appendix for completeness. 

The formulation of these boundary conditions is based upon the Riemann approach, which is illus- 
trated in figure Cl. As shown in the figure, the finite wave models assume that state 1 and state 4 are 
separated by left and right traveling acoustic waves and a convecting entropy wave. This assumption 
simplifies the mathematics for deriving the equations that were used for the boundary conditions 
because the wave structure can be solved in closed form without iteration. The finite wave boundary 
conditions have the advantage that the order of accuracy of the conditions is determined by the order of 
accuracy of the solver. 

The equations that govern the simple waves that are shown in figure Cl are (ref. 65) 


«2 - H 1 


2 n 

■y^T 


(c 2 ~ c i) 


c 

c 


2 

1 


(P 2 


Y-l 

Yy 


p i 


) 


(Cl) 


(C2) 


u 3 = «2 


(C3) 


P 3 = Pi 


(C4) 


m 4 -h 3 


2n 

Y-~l 


(C4-C3) 


f4 

C 3 


fa 


Y-l 

Yy 




(C5) 
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x 


Figure C 1 . 


Diagram of quasi-steady Riemann analysis used in derivation of finite wave boundary conditions. 
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where the subscript denotes the state, u is the flow speed, c is the sound speed, p is the pressure, y is 
ratio of specific heats, and n is a sign parameter defined to be +1 or -1 for the inflow and outflow 
boundaries, respectively (ref. 29). 

Note that equations (Cl) and (C2) describe left traveling acoustic waves, equations (C3) and (C4) 
describe convecting entropy waves, and equations (C5) and (C6) describe right traveling acoustic 
waves. 

These equations can be arranged so that states 2 and 3 are described completely by states 1 and 4: 


c 2 ~ c \ 


y — 1 

n 2~( u \ -« 4 ) + c 1 +c 4 


C,+C 4 


T~> 

fpiW 

\ P *J 


(C7) 


1 y 


Pi ~ Pi ~ P\ 


\ C \) 


Y-l 

(p^ 2 y 

y p 4j 

In 

M 3 = u 2 = u i~^zr l ( - c 2- c 0 


(C8) 


(C9) 


(CIO) 


The boundary state sought is either state 2 or state 3, depending on the direction of the flow. When 
the flow is into the domain (inflow boundary), state 2 is the correct boundary state. State 3 becomes the 
boundary state when the flow leaves the computational domain (outflow boundary). 


Example 

To illustrate the implementation of the finite wave conditions, consider the interaction of a plane 
wave with a shock with the quasi-one-dimensional Euler equations of chapter 2. 

Inflow 

The upstream exterior state is defined to be the conditions for a steady ambient flow, represented 
here by the subscript — < •». 

The acoustic wave introduced to the inflow boundary satisfies the expression 

5 p = ep ^ sin cor (C 11) 

where p_^ is the pressure at the ambient (steady) upstream state, e is an amplitude parameter, co is the 
frequency of the acoustic wave, and t is time. 

Therefore, the pressure at state 1 is 


p \ ~ P-co* 5 /* 


(Cl 2) 
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The corresponding sound speed, velocity, and density are, respectively, 


c 


l 



Y-l 

Vy 


P-o, 


(Cl 3) 


y- 1 


( C 1 — C _oo) 


(Cl 4) 


Pi = 


P l 


l 

V 


V S '-“Y 


(Cl 5) 


where j is the inflow entropy. State 4 is determined by the computational interior. Thus state 4 is 
defined by the solution in the first computational cell within the domain. 

Now that both states 1 and 4 are defined, the boundary condition is determined by solving for 
state 2 with equations (C7) through (CIO). 


Outflow 

For the outflow condition, the downstream exterior state (state 4) is defined to be steady ambient 
flow, represented by For this state, the conditions at « are determined so that the Rankine-Hugoniot 
jump conditions are satisfied across the shock, given an upstream ambient state at — The interior state 
(state 1) is defined by the computational interior, in this case the last cell of the computational domain. 
The boundary conditions are determined by solving for state 3 with equations (C7) through (CIO). 
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